i  m.  j  ria  ^  ±1  xau  - - - 

UCUWTY  CLA*Sl  F  ICATIOK  OF  THIS  fACl  f*Hmn  Data  to«^>0 


REPORT  DOCUMENTATION  PAGE 

READ  INSTRUCTIONS 

BEFORE  COMPLETING  FORM 

1.  nIpoRT  NUMBER  2.  OOVT  accuiwn  no. 

IIHR  Report  No.  304 

to  RECIPIENT*!  CATALOG  NUMBER 

4.  TITLE.  Svbtttt*) 

Viscous-Inviscid  Interaction  With  Higher-Order 
Viscous-Flow  Equations 

to  TYPE  OF  REPORT  k  PERIOD  COVERED 

Cechnical  Report 

Sept.  1983  -  Aueust  1986 

*.  PERFORNINO  ORO.  REPORT  HUMBER 

IIHR  ReDort  No.  304 

7.  AUTMORfR) 

Frederick. Stern,  Sungyul  Yoo,  and 

Virendra  C.  Patel 

to  CONTRACT  OR  GRANT  NUMSERf *) 

N0O0-14-83-K-0136 

».  PERFORMING  ORGANIZATION  NAME  AND  ADDRESS 

Iowa  Institute  of  Hydraulic  Research 

The  University  of  Iowa 

Iowa  City,  Iowa  52242 

*0.  PROGRAM  ELEMENT.  PROJECT.  TASK 

AREA  A  WORK  UNIT  NUMBERS 

NR  655-002 

It.  CONTROLLING  OFFICE  NAME  AND  ADDRESS 

Office  of  Naval  Research 

800  North  Quincy  Street 

Arlington,  Virginia  22217 

U.  REPORT  DATE 

August  1986 

IS.  NUMBER  OF  PAOC1 

Ito  MONITORING  AGENCY  NAME  A  ADDRESS^*  Mtttmrmnt  tram i  Controlling  Otttom) 

Office  of  Naval  Research 

536  South  Clark  Street 

Chicago,  Illinois  60605 

Ito  SECURITY  CLASS.  (*t  *Um  r+**X) 

Unc lassif ied 

is^  declassjfication/downoradino 

SCHEDULE 

IS.  DISTRIBUTION  STATEMENT  (of  ttit*  Rmi>ort) 

Approval  for  Public  Release;  Distribution  unlimited 

17.  DISTRIBUTION  STATEMENT  (•/  tfi*  mkmtrmet  mm MrN  M  Ml* c±  SO,  it  mttmont  trom  A^o it) 

IS.  SUPPLEMENTARY  NOTES 

IS.  KEY  WORDS  (Conttrw  on  »mM  »id«  it  m*o oooory  mnd  iMmnttty  try  WocJt  « noNrj 

Thick  3  D  Boundary  Layer,  ViscousTlnviscid  Interaction,  Partially-Parabolic 
Equations,  Computational  Fluid  Dynamics 

SO.  ABSTRACT  (Cmnttmm  m  r>rw»  mi  Mm  tt  Moxio/  mm4  IMmrtttjr  krjr  Wo ok.  mmmkof) 

The  partially-parabolic,  or  parabolised,  Navier-Stokes  equations  for 
laminar  flow,  and  the  corresponding  Reynolds  equations  for  turbulent  flow,  are 
coupled  with  an  inviscid-flow  solution  procedure  to  develop  a  viscous-inviscid 
interaction  method  which  can  be  used  in  three-dimensional  flows  which  cannot 
be  treated  by  means  of  the  classical  boundary-layer  equations.  Potential 

applications  of  such  a  higher-order  matching  procedure  are,  for  example:  thick 
boundary  layers  on  ship  stems  and  bodies  at  incidence,  interacting  shear 
layers  (wakes,  wall  jets),  solid-solid  and  solid-fluid  comers. 

DD  ,  j^Ttj  1X73  edition  or  I  NOV  •*  I*  obsolete  Unclassified 

t/n  0 20 2*0 14- 4601  |  _ 


IICUWTT  CLASSIFICATION  OF  TH»  PA0E  (Whom  D+*m  tor t«r^ 


Unclassified 

JLCU  Rl  TY  CLASSIFlCA  T1QN  OF  THIS  PAQgpTh—i  Dot*  &*tfr+d) 


..  T^'s  re.port  provides  a  detailed  overview  of  the  approach  for  general 
three-^aBMienal  flows,  and  presents  the  results  of  applications  to  some 
test  =ases.  The  Reynolds  equations  are  derived  in  nonorthogonal  curvi-- 
ni?ntr,r^+rdl^atKS-  Wlth  velocity  components  along  the  coordinate  directions, 
mP+wJ' ^+T  techm^ues*  ^ls  approach  differs  from  the  commonly-used  tensor' 
^  serves  to  establish  a  connection  with  the  more  familiar  boundary- 
layet  methods.  The  k-e  model  is  used  for  turbulent  flows.  The  partially- 
parabolic  viscous-flow  equations  are  solved  using  an  implicit  finite-differ- 
md  SIWPLER  algorithm  for  pressure-velocity  coupling.  The 
inZ^S^ld_Tl0W  solutlons  are  obtained  with  a  conforming  panel,  source-panel 
method.  Interaction  between  the  viscous  and  inviscid  regions  is  accounted  for 
using  the  displacement-body  concept.  The  relative  merits  of  interactive  and 
solutl°n  procedures  are  evaluated  by  comparing  the  viscous -inviscid 
interaction  solutions  with  large-domain  solutions  of  only  the  viscous-flow 
equations,  comparisons  are  also  made  with  experimental  data  and  other  corapu-  • 
national  methods.  Although  the  test  cases  are  restricted  to  two-dimensional 
an  axisymmetric  ilows,  the  results  clearly  demonstrate  the  feasibility  of 
higher-order  viscous-inviscid  interaction  procedures* 


Unclassified 


SECURITY  CLASSIFICATION  OF  THIS  FAOC<’Wfc*w  Dmt*  tn<«r* « 


library1 

research  reports  division 

naval  POSTGRADUATE  SCHOOL" 
MONTEREY.  CALIFORNIA  93940 


VISCOUS-INVISCID  INTERACTION  WITH 
HIGHER-ORDER  VISCOUS-FLOW  EQUATIONS, 


by 


F.  Stern,  S.  Y.  Yoo,  and  V.  C.  Patel 


Sponsored  by 

Office  of  Naval  Research 
Accelerated  Research  Initiative  (Special  Focus) 
Program  in  Ship  Hydrodynamics 

Contract  No.  N00014-83-K-0136 


IIHR  Report  304  ? 

Iowa  I  nstitute  of  Hydraulic  Research) 
The  University  of  Iowa  ? 

Iowa  City,  Iowa  52242 


August  1986 

Approved  for  Public  Release:  Distribution  Unlimited 


VISCOUS-INVISCID  INTERACTION  WITH 
HIGHER-ORDER  VISCOUS-FLOW  EQUATIONS 


by 

F.  Stem,  S.Y.  Yoo,  and  V.C.  Patel 


Sponsored  by 

Office  of  Naval  Research 

Accelerated  Research  Initiative  (Special  Focus)  Program  in  Ship  Hydrodynamics 

Contract  No.  N00014-83-K-0136 


IIHR  Report  No.  304 


Iowa  Institute  of  Hydraulic  Research 
The  University  of  Iowa 
Iowa  City,  Iowa,  52242 


August  1986 


Approved  for  Public  Release:  Distribution  Unlimited 


ABSTRACT 


The  partially -parabolic,  or  parabolised,  Navier-Stokes  equations  for 
laminar  flow,  and  the  corresponding  Reynolds  equations  for  turbulent  flow,  are 
coupled  with  an  inviscid-flow  solution  procedure  to  develop  a  viscous-inviscid 
interaction  method  which  can  be  used  in  three-dimensional  flows  which  cannot 
be  treated  by  means  of  the  classical  boundary-layer  equations.  Potential 
applications  of  such  a  higher-order  matching  procedure  are,  for  example:  thick 
boundary  layers  on  ship  stems  and  bodies  at  incidence,  interacting  shear 
layers  (wakes,  wall  jets),  solid-solid  and  solid-fluid  corners. 

This  report  provides  a  detailed  overview  of  the  approach  for  general 
three-dimensional  flows,  and  presents  the  results  of  applications  to  some 
simple  test  cases.  The  Reynolds  equations  are  derived  in  nonorthogonal  curvi¬ 
linear  coordinates,  with  velocity  components  along  the  coordinate  directions, 
using  vector  techniques.  This  approach  differs  from  the  commonly-used  tensor 
methods  but  serves  to  establish  a  connection  with  the  more  familiar  boundary- 
layer  methods.  The  k-e  model  is  used  for  turbulent  flows.  The  partially  - 
parabolic  viscous-flow  equations  are  solved  using  an  implicit  finite-differ¬ 
ence  scheme  and  the  SIMPLER  algorithm  for  pressure-velocity  coupling.  The 
inviscid-flow  solutions  are  obtained  with  a  conforming  panel,  source-panel 
method.  Interaction  between  the  viscous  and  inviscid  regions  is  accounted  for 
using  the  displacement-body  concept.  The  relative  merits  of  interactive  and 
global  solution  procedures  are  evaluated  by  comparing  the  viscous-inviscid 
interaction  solutions  with  large-domain  solutions  of  only  the  viscous-flow 
equations.  Comparisons  are  also  made  with  experimental  data  and  other  compu¬ 
tational  methods.  Although  the  test  cases  are  restricted  to  two-dimensional 
and  axisymmetric  flows,  the  results  clearly  demonstrate  the  feasibility  of 
higher-order  viscous-inviscid  interaction  procedures. 
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I.  INTRODUCTION 


The  classical  boundary-layer  equations,  which  are  based  on  the  assumption 
that  the  viscous  layer  is  thin  relative  to  the  local  radii  of  curvature  of  the 
surface,  are  parabolic  and  neglect  all  influences  of  the  downstream  flow  other 
than  those  contained  in  the  inviscid-flow  pressure  field.  Also,  in  the  clas¬ 
sical  theory,  the  inviscid  flow  is  calculated  without  accounting  for  the 
boundary-layer  displacement  effects.  In  spite  of  these  approximations,  there 
would  seem  to  be  no  question  as  to  the  tremendous  success  of  boundary— layer 
theory.  The  requirements  of  the  theory  are  met  in  many  practical  flow  situa¬ 
tions,  at  least  for  a  portion  of  the  flow  domain,  and  for  cases  where  they  are 
not  met  it  provides  a  formal  framework  upon  which  refinements  and  modifica¬ 
tions  can  be  applied  and  understood.  The  conditions  of  boundary-layer  theoiy 
are  not  met  in  a  variety  of  practical  circumstances,  for  example:  in  regions 
or  in  the  vicinity  of  flow  separation;  near  leading  and  trailing  edges;  in 
comers;  at  the  juncture  of  the  boundary  layer  and  the  free  surface  for  sur¬ 
face-piercing  bodies;  in  regions  of  strong  mass  injection;  and  in  regions  of 
strong  shock-wave  boundary-layer  interaction.  In  these  cases,  some  or  all  of 
the  terms  neglected  in  the  Navier-Stokes  equations  to  obtain  the  boundary- 
layer  equations  become  important  and,  as  a  result,  the  classical  approach 
fails  to  predict  such  flows.  There  is,  therefore,  a  need  for  approaches  which 
solve  viscous -flow  equations  which  are  more  general  than  the  boundary-layer 
equations. 

There  are  two  possible  approaches  to  the  solution  of  higher— order  vis¬ 
cous-flow  equations:  a  global  approach  in  which  one  set  of  governing  equations 
that  are  appropriate  for  both  the  inviscid-and  viscous-flow  regions  are  solved 
using  a  large  solution  domain  so  as  to  capture  the  entire  zone  of  viseous- 
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inviscid  interaction;  and  an  interactive  approach  in  which  different  sets  of 
governing  equations  are  used  for  the  viscous-  and  inviscid-flow  field  regions 
and  the  complete  solution  is  obtained  iteratively  and  interactively  through 
the  use  of  an  interaction  law,  i.e.,  patching  or  matching  conditions.  T he 
latter  approach  is  often  referred  to  as  a  zonal  approach.  It  should  be  recog¬ 
nized  that  the  former  approach  is  somewhat  more  rigorous  since  it  does  not 
rely  on  the  patching  conditions  which  are  not  exact.  Nonetheless,  for  a 
variety  of  reasons,  some  of  which  will  be  discussed  subsequently,  both  types 
of  approaches  are  of  interest.  A  complete  review  of  the  literature  is  beyond 
the  scope  of  this  report.  However,  an  overview  is  given  using  selected  refer¬ 
ences  as  examples  for  the  purpose  of  putting  the  present  work  in  perspective. 

Traditionally,  interaction  studies  have  coupled  the  thin -boundary-layer 
equations  with  inviscid-flow  solutions  which  include  viscous-flow  effects 
using  displacement -body  or  equivalent-source  methods.  Usually,  only  one 
iteration  is  performed  either  because  sufficient  accuracy  has  been  obtained  or 
due  to  slow  convergence.  Such  methods  have  been  successful  in  predicting 
flows  where  the  boundary-layer  equations  are,  in  fact,  a  good  approximation  in 
the  viscous-flow  region,  e.g.,  thin  airfoils  and  wings  at  sufficiently  small 
angles  of  attack  such  that  there  is  no  flow  separation.  Extensions  for  thick¬ 
er  airfoils  and  wings  and/or  larger  angles  of  attack,  such  that  only  a  limited 
separation  occurs  that  is  confined  to  a  thin  layer  adjacent  to  the  surface, 
have  also  been  made.  In  this  case,  the  singularity  of  the  boundary-layer 
equations  at  separation  is  removed  by  using  the  inverse  mode  and  single-pass 
solutions  can  be  obtained  using  the  FLARE  approximation.  Here  again,  the 
results  are  very  good  and  under  most  circumstances  the  interactive  boundary 
layer  procedures  can  predict  the  flow  as  well  as  the  global  Navier-Stokes  or 
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partially -parabolic  Navier-Stokes  solutions  (McDonald  and  Briley,  1983;  Mehta 
et  al.,  1985).  The  interactive  procedures  have  the  advantage  of  computational 
efficiency  over  the  global  methods. 

Applications  of  either  the  traditional  (direct-mode)  or  inverse-mode 
interactive  procedures  to  flows  in  which  the  boundary-layer  equations  do  not 
represent  a  good  approximation  (e.g. ,  axisymmetric  bodies  and  ship  hulls)  have 
generally  had  only  a  limited  success.  In  the  case  of  the  traditional  proced¬ 
ures  this  has  been  well  demonstrated  by  the  extensive  experimental  and  theor¬ 
etical  studies  of  Huang  and  associates  at  the  David  Taylor  Naval  Ship  Research 
and  Development  Center  (DTNSRDC )  (Huang  et  al. ,  1978,  1980,  1983).  For 
axisymmetric  bodies  with  relatively  sharp  trailing  edges,  the  viscous-inviscid 
interaction  is  relatively  weak,  and  the  traditional  procedures  allow  the 
boundary-layer  calculations  to  go  beyond  the  premature  separation  predicted  by 
the  classical  theory  and  show  good  agreement  with  the  experimental  data  except 
at  the  extreme  tail  region.  This  improvement  over  the  without  interaction 
solutions  is  no  doubt  due  to  the  removal  of  the  rear  stagnation  point  in  the 
inviscid-flow  solution.  However,  for  axisymmetric  bodies  with  blunt  trailing 
edges  and  more  complex  three-dimensional  bodies  the  viscous-inviscid  interac¬ 
tion  is  strong  and  the  agreement  with  the  experimental  data  has  not  been 
satisfactory.  There  have,  in  fact,  been  only  a  very  limited  number  of  in¬ 
verse-mode  interaction  studies  for  axisymmetric  and  three-dimensional  bod¬ 
ies.  Here  again,  difficulties  have  been  encountered  for  strong  interaction 
applications  (Piquet  and  Visonneau,  1985). 

Patel  (1982)  has  reviewed  the  experimental  data  for  the  viscous-flow  over 
the  stem  of  axisymmetric  bodies  and  ship  hulls,  which  he  refers  to  as  thick¬ 
boundary-layer  flows,  and  points  out  the  following  features:  (a)  for  practical 
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body  geometries  there  is  an  absence  of  flow  separation;  (b)  a  rapid  thickening 
of  the  boundary  layer;  (c)  variation  of  pressure  across  the  boundary  layer 
i nip lying  strong  viscous-inviscid  interaction;  (d)  the  development  of  a  large 
longitudinal  vorticity  component  which  may  or  may  not  lead  to  a  free-vortex 
type  separation;  (e)  a  general  reduction  in  the  level  of  turbulence.  Patel 

■  3 

concludes  that  the  appropriate  governing  equations  for  the  viscous -flow  region 
are  the  partially -parabolic  Navier-Stokes  or  Reynolds  equations.  Thus,  high¬ 
er-order  equations  must  be  used  in  the  viscous-flow  region  for  such  applica¬ 
tions. 

Investigations  using  higher-order  equations  for  axisymmetric  and  three- 
dimensional  bodies  have  been  quite  varied  in  the  approximations  embodied  and 
the  turbulence  model  and  numerical  procedures  utilized.  Generally  speaking, 
the  results  from  most  of  these  investigations  have  shown  an  improvement  over 
the  traditional  interaction  procedures;  however,  they  do  not  show  overall  good 
agreement  with  the  experimental  data.  This  can  be  attributed  to  a  variety  of 
causes:  lack  of,  inconplete,  or  incorrect  viscous-inviscid  interaction  proce¬ 
dures;  inconsistent  approximations  in  the  equations;  velocity-pressure  coupl¬ 
ing  procedures;  turbulence  modeling;  and  coordinates  and  grid  dependence.  For 
axisymmetric  bodies,  see  for  example,  Brune  et  al.  (1975),  Lee  (1978),  Dietz 
(1980),  Muraoka  (1980),  Markatos  (1984),  and  Marlin  et  al.  (1985).  For  three- 
dimensional  bodies,  see  for  example,  Abdelmeguid  et  al.  (1979),  Markatos  et 
al.  (1980),  Muraoka  (1980,  1982),  Tzabiras  and  Loukakis  (1983),  Tzabiras 
(1983),  and  Hoekstra  and  Raven  (1985).  In  most  of  the  above  references  the 
outer  boundary  was  placed  at  about  two  boundary-layer  thicknesses  from  the 
body  surface  where  conditions  are  prescribed  based  on  the  potential-flow 
solution  without  including  the  viscous-flow  displacement  effects  which  are 
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presumed  to  be  small  at  this  distance;  however,  as  pointed  out  by  Chen  and 
Patel  (1985)  and  also  by  the  present  work,  viscous-inviscid  interaction  is 
important  even  at  such  distances  from  the  body.  Thus,  these  previous  solu¬ 
tions  are  not  complete.  Also,  in  one  case  (Tzabiras  1983),  the  use  of  the 
experimental  displacement-body  to  coirpute  the  potential  flow  actually  led  to  a 
less  accurate  prediction  of  the  surface  pressure  distribution,  indicating 
other  numerical  difficulties. 

Only  Chen  and  Patel  (1985)  and  Zhou  (1982)  have  taken  a  sufficiently 
large  solution  domain  in  the  solution  of  the  partially-parabolic  Reynolds 
equations  to  capture  the  entire  zone  of  viscous-inviscid  interaction.  The  two 
methods  are  very  different.  Zhou  uses  a  streamline  iteration  method  which  may 
not  be  easily  extended  to  three-dimensional  flows.  Chen  and  Patel  use  the 
novel  finite-analytic  method  to  discretize  the  equations  and  have  presented 
results  for  both  two-  and  three-dimensional  flows.  Also,  see  the  recent 
review  article  by  Patel  and  Chen  (1985)  of  the  state-of-the-art  for  axisym- 
metric  bodies  (not  including  the  present  work).  Thus  far,  these  large  domain 
solutions  have  proven  to  be  the  most  successful. 

In  the  present  investigation,  a  viscous-inviscid  interaction  method  has 
been  developed  in  which  the  partially-parabolic,  or  parabolised,  Navier-Stokes 
equations  for  laminar  flow,  and  the  corresponding  Reynolds  equations  for 
turbulent  flow,  are  coupled  with  a  displacement-body  inviscid-flow  solution 
procedure  in  an  interactive  and  iterative  manner.  There  are  numerous  poten¬ 
tial  applications  of  such  a  higher-order  matching  procedure,  for  example: 
thick-boundary-layers  on  ship  stems  and  bodies  at  incidence,  interacting 
shear  layers,  solid-solid  and  solid-fluid  corners.  Herein,  of  particular 
interest  are  thick -boundary-layer  trailing-edge  flows,  e.g. ,  ship  boundary 
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layers.  For  such  applications  the  interaction  approach  has  the  advantage  over 
the  global  approach  in  that  it  is  most  easily  extendable  to  the  calculation  of 
ship  boundary  layers  at  nonzero  Froude  numbers. 

This  report  provides  a  detailed  overview  of  the  approach  for  general 
three-dimensional  flows,  and  presents  the  results  of  applications  to  some 
simple  test  cases.  The  Reynolds  equations  are  derived  in  nonorthogonal  curvi¬ 
linear  coordinates,  with  velocity  conponents  along  the  coordinate  directions, 
using  vector  techniques.  This  approach  differs  from  the  commonly -used  tensor 
methods  but  serves  to  establish  a  connection  with  the  more  familiar  boundary- 
layer  methods.  The  k-e  model  is  used  for  turbulent  flows.  The  partially - 
parabolic  viscous-flow  equations  are  solved  using  an  implicit  finite-differ¬ 
ence  scheme  and  the  SIMPLER  algorithm  for  pressure-velocity  coupling.  The 
inviscid-flow  solutions  are  obtained  with  a  conforming  panel,  source-panel 
method.  Interaction  between  the  viscous  and  inviscid  regions  is  accounted  for 
using  the  displacement-body  concept.  The  relative  merits  of  interactive  and 
global  solution  procedures  are  evaluated  by  comparing  the  viscous -inviscid 
interaction  solutions  with  large-domain  solutions  of  only  the  viscous-flow 
equations.  Comparisons  are  also  made  with  experimental  data  and  other  compu¬ 
tational  methods.  Although  the  test  cases  are  restricted  to  two-dimensional 
and  axisymmetric  flows,  the  results  clearly  demonstrate  the  feasibility  of 
higher-order  viscous-inviscid  interaction  procedures. 


II,  METHOD  OF  APPROACH 

Consider  the  flow  field  in  the  vicinity  of  a  body  fixed  in  a  uniform 

stream  with  velocity  U  of  an  incompressible  viscous  fluid.  The  body  shape  is 

o 
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assumed  to  be  sufficiently  streamlined  and  the  Reynolds  number  sufficiently 
large  such  that  no  flow  separation  occurs  and  viscous  effects  are  confined  to 
a  relatively  narrow  region  adjacent  to  the  body  surface  and  in  the  wake.  As 
depicted  in  figure  1,  the  flow  field  can  be  divided  into  three  regions. 
Region  1  is  the  inviscid-flow  region.  Region  2  is  the  thin-boundary'-layer 

which  ends  at  a  station  xu  beyond  which  boundary -layer  approximations  are  no 
longer  valid.  Region  3  (X  >  xu)  is  the  thick-boundary-layer  and  wake  region 
in  which  it  is  assumed  that  only  streamwise  gradients  of  viscous  and  turbulent 
stresses  can  be  neglected.  The  inviscid-flow  region  extends  to  a  distance  yQ 
beyond  which  uniform  stream  conditions  may  be  assumed.  Appropriate  computa¬ 
tional  methods  can  be  used  for  each  of  the  three  regions.  They  are  related 
through  their  boundary  conditions  and  are  not  necessarily  independent.  In 
particular,  the  flow  fields  in  regions  3  and  1  are  interdependent  such  that  a 
complete  solution  for  each  region  has  to  be  determined  iteratively  through  the 
use  of  a  viscous-inviscid  interaction  procedure.  Alternatively,  the  flow 
fields  in  regions  3  and  1  can  be  solved  simultaneously  by  simply  extending 
region  3  to  also  include  the  portion  of  region  1  influenced  by  the  interac¬ 
tion.  Such  a  large-domain  solution  captures  the  entire  zone  of  viscous- 
inviscid  interaction  and  can  be  used  to  assess  the  accuracy  of  a  small-domain 
interaction  solution.  Herein,  both  types  of  solutions  will  be  obtained  in 
order  to  explicate  the  nature  of  an  interactive  solution. 


III.  VISCOPS-INVISCID  INTERACTION 

In  an  interactive  approach  to  the  present  problem,  the  viscous-  and 
inviscid-flow  regions  3  and  1  are  demarcated  by  the  common  boundary  6  ,  as 
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shown  in  figure  1.  This  boundary  must  be  placed  at  a  sufficient  distance  from 
the  body  surface  and  wake  centerplane  such  that,  exterior  to  6  (region  1), 
viscous  effects  can  be  neglected.  Traditionally,  this  boundary  is  placed  at 
or  just  beyond  the  edge  of  the  boundary  layer.  Note  that  this  should  be 
considered  as  the  minimum  value  for  6  ,  and  it  is  also  possible  to  place  <5  at 
distances  greater  than  the  boundary -layer  thickness,  as  has  been  done  by  other 
investigators.  Herein,  the  traditional  definition  is  used.  The  flow-field 
solutions  in  regions  3  and  1  are  obtained  separately  and  by  quite  different 
means,  but  they  are  interdependent  through  the  common  boundary  condition  that 
the  two  solutions  merge  smoothly  at  6 .  That  is,  the  viscous-flow  solution 
(region  3)  is  obtained  with  edge  conditions  at  6  specified  based  on  the  invis- 
cid-flow  solution  (region  1)  which  is  obtained  with  recognition  of  the  dis¬ 
placement  effect  of  the  viscous  flow.  The  complete  solution  is  obtained 
iteratively  until  convergence  is  achieved. 

It  has  long  been  recognized  that  the  viscous-flow  region  displaces  the 
inviscid-f low  streamlines  such  that  the  inviscid  flow  is  not  the  same  as  that 
about  the  actual  body,  but  rather  that  about  a  surface  displaced  into  the 
fluid  a  distance  6  referred  to  as  the  displacement  surface.  The  displacement 
surface  can  be  defined  unambiguously  by  the  following  two  requirements:  (a) 
that  it  be  a  stream  surface  of  the  inviscid  flow  continued  from  outside  the 
boundary  layer;  (b)  that  the  inviscid-f low  discharge  between  this  surface  and 
any  stream  surface  exterior  to  the  boundary  layer  be  equal  to  the  actual 
discharge  between  the  body  and  the  latter  stream  surface.  The  second  condi¬ 
tion  implies  that  the  flow  reduction  inside  the  viscous  flow  is  compensated  by 
an  outward  displacement  of  such  a  stream  surface  through  a  distance  6  ,  i.e. 
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(III-l) 


/  v  .  dA  =  /  (V  -  V)  .  dA 

A  — £  A  _E 

AS*  A6 

where  is  the  velocity  vector  of  the  outer  inviscid  flow  analytically  con¬ 
tinued  into  the  viscous-flow  region,  V  is  the  viscous-flow  velocity  vector, 
and  Ag*  and  are  the  cross-sectional  areas  between  the  actual-body  surface 
and  the  displacement-body  surface  and  the  boundary-layer  surface  respec¬ 
tively.  Thus,  the  inviscid— flow  solution  is  obtained  for  the  displacement 
body .  This  solution  then  provides  the  boundary  conditions  for  the  viscous 
flow 

U(6  )  =  U  (5  )  =  U 
P  e 

W(6  )  =  W  (5  )  =  W  ( III-2) 

p  e  ' 

p(<S )  =  p  (6  )  =  p 
P  e 

Since  6  and  V  (6  )  are  not  known  a  priori,  an  initial  guess  must  be  provided 
and  the  complete  solution  obtained  by  iteratively  updating  the  viscous-  and 
inviscid-flow  solutions  until  the  "patching"  conditions  (III-l)  and  (IH-2) 
are  satisfied.  Note  that,  in  the  present  viscous-inviscid  interaction  pro¬ 
cedure,  no  assumptions  have  been  made  with  regard  to  the  thickness  of  the 
boundary  layer.  The  primary  assumption  is  that  the  inviscid—  and  viscous— flow 
solutions  can  be  patched  together  through  conditions  (III-l)  and  (III-2). 

Within  the  context  of  thin-boundary-layer  theory,  and  for  two-dimensional 
flow,  equation  (III-l)  becomes 
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( I II—  3 ) 


where 


«*  =  -jj-/  (1  "  U/Ue}  ^ 

e  o 

U  =  U  (o)  =  U  (<$  ). 
e  p  p 


Lighthill  (1958)  refers  to  this  definition  of  S*  (III-3)  as  the  flow-reduction 
method  which  he  shows,  for  both  two-  and  (with  the  appropriate  definition) 
three-dimensional  flow,  is  completely  equivalent,  within  the  context  of  thin- 
boundary-layer  theory,  to  three  other  definitions:  equivalent  source,  velo¬ 
city  comparison,  and  mean  vorticity.  The  velocity-comparison  method  was  first 
introduced  by  Moore  (1953)  and  is  closely  related  to  the  equivalent -source 
method,  in  which  it  is  shown  that  the  displacement  effect  of  the  boundary 
layer  on  the  inviscid  flow  can  be  represented  by  an  additional  distribution  of 
sources  on  the  actual  body  surface  of  strength 


a 


BL 


1_ 

2ir 


d_ 

dx 


:u/> 


( III-4) 


where  x  is  the  distance  along  the  body  surface  in  the  streamwise  direction. 

Landweber  (1978)  has  pointed  out  that  ( III— 4)  is  just  the  first  approximation 

to  the  solution  of  the  general  integral  equation  for  a_.T  due  to  the  boundary- 

d.  ^ 

layer  outflow  velocity  distribution  -r—  (Ue6  ).  To  the  same  order  of  approxima¬ 
tion,  the  source  distribution  for  the  actual  body  is  given  by 

hs-irhW  (m-5) 

where  n^  is  the  x-eomponent  of  the  unit  normal  to  the  body  surface.  Based  on 
either  thin-  or  thick-boundary-layer  order-of -magnitude  estimates  (see  Table 
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1)  cr_T  /  U  ~  0  (e )  where  e  is  the  nondimensional  boundary-layer  thickness 
BL  o 

and  aB  /U  ~  0(n1).  The  equivalent-source  method  has  been  used  extensively 
and  with  success  in  viscous-inviscid  interaction  procedures  for  airfoils  and 
wings  where  n^  is  small  in  the  region  of  interaction.  However,  for  the  prob¬ 
lems  of  present  interest  (e.g.  axisymmetric  bodies  and  ship  hulls),  n^  is  not 
necessarily  small  near  the  trailing  edge  and  the  equivalent-source  method  is 
not  useful  for  representing  the  viscous-flow  displacement  effect  on  the  invis- 
cid  flow. 


IV.  VISCOUS  FLOW 

In  the  thick  boundary  layer  and  wake  (region  3  in  figure  1)  it  is  assumed 
that  only  streamwise  gradients  of  viscous  and  turbulent  stresses  can  be  ne¬ 
glected.  Under  this  assumption,  the  Reynolds  equations  are  reduced  to  a 
simplified  form  referred  to  as  the  partially -parabolic  (Reynolds)  equations. 
In  these  equations,  the  velocity  field  is  elliptic  in  transverse  planes  and 
parabolic  in  the  streamwise  direction  while  the  pressure  field  is  fully  ellip¬ 
tic.  Solutions  to  the  partially -parabolic  equations  can  be  obtained  itera¬ 
tively  by  solving  the  parabolic  equations  that  result  when  the  pressure  field 
is  specified  and  subsequently  updating  the  pressure  field  using  the  results 
from  the  latest  parabolic  solution.  Of  crucial  importance  is  the  manner  in 
which  the  velocity  and  pressure  fields  are  coupled.  Many  procedures  have  been 
tried  for  this  purpose  (e.g.,  Anderson  et  al.  1984).  In  the  present  work,  a 
modified  form  of  the  SIMPLER  algorithm  (Patanker,  1980)  that  enhances  global 
convergence  has  been  used  (Chen  and  Patel,  1985).  Selection  of  the  appropri¬ 
ate  coordinate  system  and  grid  generation  technique  used  to  obtain  the  coor- 
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dinates  is  also  important.  A  streamline  coordinate  system  is  the  most  consis¬ 
tent  with  the  assumptions  of  the  partially-parabolic  equations;  however,  such 
a  coordinate  system  is  difficult  to  generate.  Alternatively,  body-fitted 
coordinates  can  be  used  in  which  the  axial  coordinate  should  be  roughly 
aligned  with  the  streamlines  since  they  are  coincident  on  the  body  surface 
itself.  The  partially-parabolic  assumptions  are  made  in  this  preselected 
axial  coordinate  direction.  In  the  present  work,  both  simplified  analytic  and 
numerically  generated  body-fitted  coordinate  systems  have  been  used.  The 
Reynolds  stresses  are  modeled  using  the  k-e  turbulence  model.  A  complete 
transformation  of  the  governing  equations  is  used  such  that  the  directions  of 
the  velocity  components  are  along  the  grid  lines.  It  should  be  recognized 
that  very  few  investigators  have  used  a  complete  transformation  of  the  govern¬ 
ing  equations,  no  doubt  due  to  the  complexity  of  their  derivation  as  will 
become  apparent  subsequently.  The  more  common  approach  is  a  partial  transfor¬ 
mation  in  which  only  the  coordinates  are  transformed  and  the  velocity  compon¬ 
ents  are  maintained  in  either  cartesian  or  polar  coordinates.  The  governing 
equations  are  reduced  to  algebraic  form  using  finite  differences  and  solved 
implicitly  by  the  method  of  lines.  In  the  subsections  to  follow,  the  details 
of  these  computational  procedures  are  discussed. 

A.  Equations  and  Coordinate  ^ysteaa.  The  partially-parabolic  equations 
are  solved  using  a  nonorthogonal  curvilinear  coordinate  system  in  which  the  x- 
coordinate  is  roughly  aligned  with  the  flow  direction  and  the  y-coordinate  is 
in  a  plane  transverse  to  the  body  axis  X  (see  figure  la).  For  three-dimen¬ 
sional  flow,  the  z-coordinate  is  also  in  the  transverse  plane  and  in  the 
girthwise  direction  (see  figure  lb).  The  Reynolds  equations  in  nonorthogonal 
curvilinear  coordinates  can  be  derived  either  through  the  use  of  vector  or 
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tensor  analysis.  Here,  a  vector  approach  has  been  used  since  it  lends  itself 
to  more  physical  insight.  Recently,  Richmond  et  al.  (1986)  have  provided  a 
derivation  using  tensor  analysis. 

The  Reynolds,  continuity,  turbulent-kinetic-energy  k,  and  its  dissipa¬ 
tion-rate  e  equations  for  steady  incompressible  flow  can  be  written  in  the 
following  vector  form: 


2  V  (V*  V)  -  V  x  co_  =  -  V  p/p  +  v{V  (V*  V)  -  Vxw} 

-  V  •  v,v.  +  (v. )  V»  v 
i  J  i  - 


(IV-1) 


V*V  =  0 


(IV-2: 


V*Vk  =  V-  (—  Vk)  +  G  -  e 
ak 


V-7C  =7  •  (^7e)  *  Cel  G  |  -  0e2  |2 


(IV-3) 


(IV-4) 


where _V  =  (U,V,W)  are  the  mean  velocity  components,  _v  =  (u,v,w)  are  the  turbu¬ 
lent  velocity  components,  p  is  the  mean  pressure,  =  V  x  V  is  the  mean  vorti- 


C-fty,  Vivj  ^re  ^®y^dds  stresses  (the  overbar  denotes  time  averaging), 


k  =  -k  v* v  is  the  turbulent  kinetic  energy,  v.  =  C  k  /e  is  the  eddy  visco- 
sity,  and  G  is  the  turbulence  generation.  Since  the  fluid  is  assumed  to  be 
incompressible,  the  terms  involving  V  •  V  and  V  •  v  in  equation  (IV-1)  are 
identically  zero,  but  have  been  included  since  they  aid  in  putting  the  trans¬ 
formed  equations  into  a  more  compact  form.  The  usual  values  are  used  for  the 
constants  in  the  k-e  equations,  namely,  (C  ,  a  a  ,  C  C  )  =  (.09  1. 

1.3,  1.44,  1.92).  The  turbulence  generation  term  is  defined  by 
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( IV- 5) 


G  vt[2(En  +  e22  +  e^)  +  4  (e12  +  £33  +  e3i}1 
where  e  is  the  rate -of -strain  tensor 

eij  =  \  IV V  +  VVTJ  ( IV— 6) 

T 

In  (IV-6)  VV  is  the  deformation-rate  tensor  and  VV  its  transpose,  i.e. 

T 

VV  =  .  The  Reynolds  stresses  required  in  (IV-1)  are  related  to  k  and  e 

through  the  isotropic  eddy  viscosity  concept: 


V  V 

ij 


=  -  2ViJ 


+  3  k(hihjglj 


( IV-7) 


where  the  h^  are  the  metric  coefficients  and  is  the  inverse  metric  tensor 
both  of  which  are  defined  below. 

Equations  (IV-1)  -  (IV-7)  can  be  transformed  into  any  coordinate  system 
through  the  use  of  appropriate  definitions  of  the  gradient  (V ) ,  divergence 
(V*  ),  and  curl  (Vx)  vector  operators.  The  details  of  this  procedure  for 
orthogonal  curvilinear  coordinates  are  provided  by  Rouse  (1959).  For  nonorth- 
ogonal  curvilinear  coordinates  the  appropriate  vector  operator  definitions  are 
not  readily  available.  They  were  probably  first  derived  by  Weatherbum 
(1926).  Following  Weatherburn,  and  referring  to  figures  1  and  2  for  the  pre- 

A  A  A  A 

sent  notation,  the  unit  vectors  e^^  =  (e1,e2,e^)  in  the  directions  of  the 
nonorthogonal  curvilinear  coordinates  (x,y,z)  are  defined  in  terms  of  the  body 
cartesian-coordinate  position  vector 
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( IV— 8 ) 


A  A 

R  =  X  (x,y,z)  i  +  Y  (x,y,z)  j  +  Z  (x,y,z)  k 


by 


6-j^  Rx/hp  fiy/hg >  Rz/h^  ( IV-9) 

where 

hl  =  \Rx\>  h2  =  lRy  I  >  h3  =  lRJ  CIV-10) 

and  a  lettered  subscript  denotes  a  partial  derivative.  The  angles  (A  ,p  ,v  ) 
between  the  (x,y,z)  coordinate  axes  are  given  by 

A  A 

cos  X  =  e^* 

A  A 

cos  p  =  e1*  e^  (IV-11) 

cos  v  =  e^* 


and  the  unit  normals  to  constant  x-  y-  and  z-surfaces  are  given  respectively 

by 


A  A 


^2  ^  e3  —  lih  s  ^^2^2  ^  Gh^e^  ] 

2  3 


A 


e3  X  el  =  h^F  [Hhlel  +  Bh2e2  +  Fh3e3I 


( IV-12) 


el  X  e2  Il^s  [Ghlel  +  Fh2e2  Ch3e3] 
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where  s  is  the  triple  product 


A  A 


s  =  (h^l^h^)  (e^*  e2  x  e^) 


[Ah-^  +  Hh^t^  cos  v  +  Gh^h^  cos  W  1 


1/2 


( IV— 13 ) 


and 


.  ,Z2  .  2, 

A  =  t^h^  sin  X 

2  2  2 

B  =  h-^h^  sin  u 

2  2  2 

C  =  h^i^  sin  v 

F  =  (h^h^  cos  m )  (h.^  cos  v) 

G  =  (h-jh2  cos  v )  (h2h^  cos  X ) 

H  =  (h2h^  cos  X)  (hjh^  cos  jj  )  -  h^Oi^t^  cos  v) 


h^  (t^h^  cos  X  ) 
2 

(h^h^  cos  u ) 


( IV— 14 ) 


The  inverse  metric  tensor  is  defined  by 


“WYj1 


-1 


A  H  G 
7  [  HBF 
G  F  C 


(IV-14.1) 


In  terms  of  the  above  quantities,  the  gradient  of  any  scalar  Q(x,y,z)  and 

AAA 

divergence  and  curl  for  any  vector  V(x,y,z)  =  +  Ve  are  given  by: 


VQ  =  —  {  (AQx  +  HQy  +  GQz)  h ^  (HQ^.  +  BQ^  +  FQz ) 
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(IV-15) 


+  (GQx  +  FQy  +  CQz)  h3e3} 


1  3  .  sV_  a  sV„ 

V*V  =  —  f—  ( - -)  +  —  ( _ -)  +  i_  I 1 

-  s  3 x  [h1  9y  lh2  '  9z  lh 


)] 


(IV-16) 


1  r  9 


VxV  -  —  {gy  [V-^h3  cos  y  +  V2h3  cos  X  + 


3z  ^1^2  cos  v  +  ^2^2  +  V3h2  cos  ^lel 


1  r8 


+  "s  ^az"  ^1^1 +  ^2^1 cos  v  +  ^h^cos  M 


-  9^  fVih3C0S  w  +  Vghysos  X  +  vy^]}  h2e2 


1  3 

+  "s  ^a7  [vih2cos  v  +  V2h2  +  V3h2  cos  X 


“  9^  ^Vlhi+  V2hlcos  v  +  V3hicos  U ]}  h3e3 


( IV-17) 


The  transformed  equations  are  very  lengthy  and  are  provided  in  Appendix 
I.  The  equations  have  been  put  in  a  form  that  is  similar  to  that  used  by  Nash 
and  Patel  (1972)  for  orthogonal  curvilinear  coordinates.  By  comparison,  it  is 
seen  that,  for  the  present  circumstances,  the  coefficients  in  the  governing 
equations  depend  on  terms  related  to  not  only  the  curvatures  of  the  coordin¬ 
ates  but  also  their  angular  orientation.  Due  to  the  complexity  of  the  deriva¬ 
tion  of  the  transformed  equations  it  was  desired  to  check  their  accuracy; 
however,  this  was  made  difficult  by  the  fact  that  no  other  presentations  of 
the  governing  equations  in  nonorthogonal  curvilinear  coordinates  in  a  format 
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and  notation  similar  to  the  present  one  are  known  to  exist.  The  following 
checks  were  made:  for  (A,u,v)  90°,  the  orthogonal  form  of  the  equations  was 
recovered;  for  (A  ,v )  ->•  90°,  and  subject  to  the  boundary-layer  assumptions  the 
boundary-layer  equations  of  Cebeci  et  al.  (1978)  were  recovered;  and  some  of 
the  coefficients  were  compared  with  their  corresponding  counterparts  in  the 
tensor  form  of  the  equations  presented  by  Richmond  et  al.  (1986). 

The  partially -parabolic  equations  are  obtained  from  the  complete  set  of 
equations  in  Appendix  I  through  the  use  of  the  order-of-magnitude  estimates 
given  in  Table  1. 


Table  1:  Order-of -Magnitude  Estimates  for  a  Thick  Boundary  Layer 


Quantity 


Order-of -Magnitude 


U 

V 

W 

3/9x 

3/3y 

3/3z 

v 


1 

e 

e 

1 

e  -1 

e-1 

2 

£ 


v.  v .  e 

i  J 


Patel  (1982)  has  shown  that  these  order-of-magnitude  estimates  are  consistent 
with  the  partially -parabolic  assumptions.  Also,  in  obtaining  the  partially- 
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parabolic  equations  no  assumptions  are  made  with  regard  to  geometrical  quanti¬ 
ties;  that  is,  all  geometrical  quantities  are  considered  0(1).  The  x-momentum 
equation  is  obtained  by  retaining  terms  of  0(1)  only.  All  other  equations  are 
obtained  by  retaining  terms  of  0(1)  +  0(e).  The  partially -parabolic  equations 
are  provided  in  Appendix  II.  Lastly,  it  should  be  mentioned  that,  strictly 
speaking,  the  complete  equations  can  also  be  rendered  partially-parabolic  by 
simply  only  neglecting  the  viscous-  and  turbulent-diffusion  terms  (second- 
order  derivatives)  in  the  x  direction. 

B.  Discretization.  The  governing  equations  are  reduced  to  algebraic 
form  by  approximating  all  the  spatial  derivatives  by  finite  differences.  A 
s^&6£®red-grid  system  has  been  used  in  order  to  avoid  difficulties  in  the 
velocity-pressure  coupling  procedure  to  be  discussed  subsequently.  The  grid 
arrangement  and  notation  are  as  shown  in  figure  3*  An  implicit  finite-differ¬ 
ence  scheme  is  used  which  is  basically  only  first-order  accurate;  however, 
certain  derivatives  have  been  evaluated  using  second-order  central  differences 
and  all  terms  have  been  evaluated  at  the  proper  grid  location  by  using  aver¬ 
ages  of  the  surrounding  values.  The  detailed  finite-difference  formulas  are 
provided  in  Appendix  III.  The  overall  procedure  is  described  below. 

The  x— momentum  equation  is  discretized  with  reference  to  control  volume  I 

of  figure  3.  Referring  to  the  partially-parabolic  form  of  the  x-momentum 

equation  (B-l),  figure  3,  and  Appendix  III. A:  only  the  x-  and  y-derivatives 

in  the  convective  acceleration  and  the  y-diffusion  terms  are  maintained  on  the 

left  hand  side  to  form  the  tridiagonal  matrix;  the  remaining  terms  are  grouped 

into  a  single  source  term  and  are  put  on  the  right  hand  side;  and  all  terms 

are  evaluated  at  the  location  of  if' 

m,n 
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The  y-momentum  equation  is  discretized  with  reference  to  control  volume 

II  of  figure  3.  Referring  to  the  partially-parabolic  form  of  the  y-momentum 
equation  (B-2),  figure  3,  and  Appendix  III. A:  only  the  x-  and  y-derivatives  in 
the  convective  acceleration  and  y-diffusion  terms  are  maintained  on  the  left 
hand  side  to  form  the  tridiagonal  matrix;  the  remaining  terms  are  grouped  into 
a  single  source  term  and  are  put  on  the  right  hand  side;  and  all  terms  are 
evaluated  at  the  location  of  v 

m,n 

The  z-momentum  equation  is  discretized  with  reference  to  control  volume 

III  of  figure  3.  Referring  to  the  partially-parabolic  form  of  the  z-momentum 

equation  ( B— 3 ) ,  figure  3,  and  Appendix  III. A:  only  the  x-  and  z-derivatives 

in  the  convective  acceleration  and  z-diffusion  terms  are  maintained  on  the 

left  hand  side  to  form  the  tridiagonal  matrix;  the  remaining  terms  are  grouped 

into  a  single  source  term  and  are  put  on  the  right  hand  side;  and  all  terms 

are  evaluated  at  the  location  of  v/ 

m,n 

The  k-e  equations  are  discretized  with  reference  to  control  volume  IV  of 
figure  3.  Referring  to  the  partially-parabolic  form  of  the  k-e  equations  (B- 
4)  and  (B-5)  respectively,  figure  3,  and  Appendix  III.B:  only  the  x-  and  y- 
derivatives  in  the  convective  acceleration  and  y-diffusion  terms  are  main¬ 
tained  on  the  left  hand  side  to  form  the  tridiagonal  matrix;  the  remaining 
terms  are  grouped  into  a  single  source  term  and  are  put  on  the  right  hand 

a  £ 

side;  and  all  terms  are  evaluated  at  the  location  of  k  and  e 

m,n  m,n 

By  means  of  the  above  finite  difference  scheme,  the  three  momentum  and 
the  k-e  turbulence-model  equations  can  be  put  in  the  form: 


a/  ,  + 

1  m-l,n 


a^l/'  + 

2  m,n 


£L.y  , 

3  m+l,n 


Sn  +  Pu 
m,n  m,n 


(IV-18) 
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(IV-19) 


( IV-20) 


( IV-21) 


m.n 


(IV-22) 


where  the  S{>  terms  are  the  source  terms  and  Fty  are  the  pressure-gradient 
terms  (if>  -  (U, V,W,k,e ) ) .  Note  that  equations  (IV-18)  -  (IV-22)  are  nonlinear 
since  both  the  coefficients  a^  through  e^  and  the  source  terms  are  functions 
of  the  unknowns  (U,V,W,k,e )  (see  Appendix  III).  If  the  pressure  field  is 
known  then  equations  (IV-18)  -  (IV-22)  can  be  solved  directly  for  the  velocity 
field  (U,V,W)  and  turbulence  parameters  (k,e)j  however,  since  the  pressure 
field  is  unknown,  it  must  be  determined  such  that  the  continuity  equation  is 
also  satisfied.  The  velocity-pressure  coupling  procedure  is  the  subject  of 
the  next  section. 

C.  Velocity-Pressure  Coupling.  The  coupling  of  the  velocity  and  pres¬ 
sure  fields  is  accomplished  through  the  use  of  a  two-step  iterative  procedure 
involving  the  continuity  equation.  In  the  first  step,  the  solution  to  the  mo¬ 
mentum  equations  for  a  guessed  pressure  field  is  corrected  at  each  cross  sec¬ 
tion  such  that  continuity  is  satisfied.  However,  in  general,  the  corrected 
velocities  are  no  longer  a  consistent  solution  to  the  momentum  equations  for 
the  guessed  p.  Thus,  the  pressure  field  must  also  be  corrected.  In  the  sec¬ 
ond  step,  the  pressure  field  is  updated  again  through  the  use  of  the  continu¬ 
ity  equation.  This  is  done  after  a  complete  solution  to  the  velocity  field 
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has  been  obtained  for  all  cross-sections  in  a  manner  that  properly  accounts 
for  the  elliptic  nature  of  the  pressure  field  and  also  accelerates  conver¬ 
gence.  Repeated  global  iterations  are  thus  required  in  order  to  obtain  the 
converged  solution.  The  details  of  the  derivation  of  the  pressure-correction 
and  pressure  equations  are  provided  in  Appendix  III.  C.  The  overall  procedure 
is  described  below. 

Both  the  pressure-correction  and  pressure  equations  have  the  same  form 
and  are  derived  in  the  same  manner  from  the  discretized  form  of  the  continuity 
equation,  i.e. 


<vL,  vi  _) 


— - -  (u^  _  7j^  i)  + - 1 -  _  v 

,  x£  m.n  m.n  .£-1/2  m+l,n  m,n 

(Ax-)  ^  (Ay-)  1/9 

m,n  m+l/2,n 


-  o 


(Az-) 


Si- 1/2  m,n+l  m,n 

m,n+l/2 


(IV-26) 


where 


U2  =  (|_)2  u2 

m,n  h^  m,n  m 

V2  =  (^-)2  / 

m,n  h^  m,n  m 

W2  =  (^-)2  / 


( IV-27) 


,n’  m,n’ 

from  the  discretized  form  of  the  momentum  equations  (IV-l8)-(IV-20)  by  putting 
equations  ( IV-18) - ( IV-20)  in  the  form 


The  velocity  components 


(if 


m 


W2  )  required  in  (IV-27)  are  obtained 
m.n 
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u*  =  u* 


(IV-28) 


where  C  n  L  /  are  referred  to  as  pseudovelocities. 

Jll  y  XI  y  III  y  II  «  III  |  Xl 


Substituting  (IV-28)  and  (IV-27)  into  (IV-26)  and  representing  the  pres¬ 
sure  gradient  terms  in  (IV-28)  by  finite  differences  results  in  the  desired 


equation  for  pressure.  In  the  IT  n  equation,  a  forward  difference  is  used  for 

p  and  central  differences  for  p.r  and  p„.  In  the  v  equation,  a  backward 
A  y  *  m,n 

difference  is  used  for  p  and  central  differences  for  pY  and  p„.  In  the  / 

y  x  z  nun 


equation,  a  backward  difference  is  used  for  pz  and  central  differences  for  px 
and  py  All  terms  are  evaluated  at  the  respective  velocity  component  location 
by  using  averages  of  the  surrounding  values.  The  resulting  equation  for 
pressure  involves  20  of  the  27  nodes  corresponding  to  x  indices  U-l,  £, 
t+1)  and  y  and  z  indices,  (m-1,  m,  m+1)  and  (n-1,  n,  n+1),  respectively ; 
however,  only  certain  terms  are  maintained  on  the  left  hand  side  to  subse¬ 
quently  form  the  tridiagonal  matrix,  and  all  remaining  terms  are  grouped  into 
a  single  source  term  and  put  on  the  right  hand  side: 


( IV-29) 
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It  should  be  recognized  that  no  approximations  have  been  made  in  deriving  (IV- 
29)  and  as  such  it  is  an  exact  representation  of  the  equation  of  continuity  in 
terms  of  pressure. 

Equation  (IV-29)  is  first  used  to  correct  the  velocity  field  obtained 
from  the  solution  of  the  momentum  equations  for  a  guessed  pressure  field.  In 

A  A 

this  case,  p  is  designated  as  p  and  both  upstream  and  downstream  values  of  p 
are  neglected.  The  resulting  pressure -correction  equation  is 


foPm  „ + 


2pm,n  +  4^ni+l,n  +  ^f^m-ljn  +  *6^111^+1 


+  f. 


+  frrP  T  =  Sp 
7^m,n-l  * 


(IV- 30) 


where  Sp  is  obtained  by  evaluating  Sp  (C-61)  with  the  current  solution  to  the 

momentum  equations  U*  ,  /  .  /  substituted  for  t  .  f  .t  .  -me 

m,n’  m,n’  m,n  m,n’  m,n’  m,n 

pressure-correction  equation  (IV-30)  is  approximate  since  both  upstream  and 

A 

downstream  values  of  p  have  been  neglected  as  well  as  the  influence  of  pres¬ 
sure  corrections  on  the  neighboring  velocities.  The  latter  approximation 
neglects  the  indirect  or  implicit  influence  of  the  pressure  correction  on 
velocity  and  is  the  reason  for  the  words  semi-implicit  in  the  name  SIMPLER. 
Note  that  both  approximations  are  justified  since  the  corrected  velocity  field 

A 

is  only  an  intermediate  solution,  and  when  the  solution  converges,  p  is  zero. 
Equation  (IV-30)  is  solved  at  each  cross-section  l  using  the  method  of  lines 
for  two-dimensional  and  axisymmetric  flow  applications  and  an  alternating 
direction  implicit  method  (ADI)  for  three-dimensional  flow  applications. 

A 

Subsequently,  the  velocity  field  is  corrected  using  equations  (IV-28)  with  p 

substituted  for  p  and  with  the  current  solution  to  the  momentum  equations 

if'  ,  V^  ,l/  substituted  for  l/  V^  l/’ 
m,n’  m,n’  m,n  m,n>  m,n’  m,n 
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When  the  exit  cross-section  £=LL  is  reached,  equation  (IV-29)  is  again 
used  to  update  the  pressure  field.  In  this  case,  no  approximations  are  made. 
The  pressure  equation  is  solved  by  marching  from  downstream  to  upstream  using 
the  method  of  lines  for  two-dimensional  and  axisymmetric  flow  applications  and 
an  ADI  method  for  three-dimensional  flow  applications  at  each  cross-section  £ 
with  the  £+1  and  £-1  values  considered  as  known  from  the  previous  iteration. 
With  a  new  pressure  field  thus  obtained,  the  entire  process  is  repeated  until 
the  results  converge;  that  is,  a  compatible  velocity-  and  pressure-field 
solution  is  found. 

.2; _ Boundary  Conditions .  The  partially -parabolic  equations  require 

boundary  conditions  for  the  pressure,  either  explicit  or  implicit,  on  all 
boundaries.  Boundary  conditions  are  required  only  from  upstream  and  in  the 
cross-section  for  the  velocity  coirponents  and  turbulence  parameters.  Note 
that  only  three  of  the  unknowns  (U,V,W,p)  can  be  specified  on  each  boundaxy^or 
the  problem  is  over  specified.  The  fourth  unknown  is  determined  through  the 
solution  of  the  governing  equations. 

Referring  to  figure  4,  the  boundary  conditions  used  in  solving  the  momen¬ 
tum  and  k-e  equations  (IV-l8)-(IV-22),  and  the  pressure  (IV-29)  and  pressure- 
correction  (IV-30)  equations  are  as  follows: 

Inlet  Boundary  Sj 

Initial  conditions  for  (U,V,W,k,e)  are  specified  based  on  simple  flat- 

A 

plate  solutions.  Initial  conditions  for  p  and  p  are  not  required. 

Exit  Boundary  Sg 
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A  zero-gradient  condition  is  used  for  p  and  p  and  no  conditions  are 


required  for  the  velocity  components  and  k-e  . 


Outer  Boundary  Sp 


U  =  U  ,  W  =  W  ,  p=p 
e'  e’  r  *e 


(IV-31) 


where  (U^,  Wg,  pg)  are  specified  from  the  inviscid-flow  solution.  Large- 
domain  solutions  are  obtained  by  placing  the  outer  boundary  at  a  sufficient 
distance  from  the  body  surface  and  simply  specifying  (Ue,  We,  pe)  =  (1,  0., 

0.  ). 

Body  Surface  Boundary  Sp 

For  laminar  flow,  the  solution  is  carried  out  up  to  the  body  surface 
where  the  no-slip  condition  is  applied:  U  =  V  =  W  =  0.  For  turbulent  flow, 
the  wall-function  approach  of  Chen  and  Patel  (1985)  is  used.  In  this  proced¬ 
ure,  the  first  two  grid  points  are  placed  in  the  log-law  region.  With  a 
guessed  value  of  the  wall-shear  velocity  U  the  required  boundary  conditions 
at  the  first  grid  point  for  the  velocity  components  (U,V,W)  and  k  and  e  are 
obtained  from  the  log-law  and  the  assumption  of  local  equilibrium: 


]  +  2[  (1  +  ATy  +  )1^2  -  1]} 


+  B  +  3.7  A 


(IV-32) 


P 
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( IV— 33 ) 


k  =  U//c 

T  U 


e  =  UT/<y 


where  U  is  the  wall-shear  velocity  defined  by  U  =  /  t  /p ,  y+  =  yU  /v  is  the 


w 


dimensionless  distance  measured  in  the  direction  normal  to  the  surface, 

A?  =  is  the  dimensionless  pressure  gradient,  A  is  the  dimensionless 

P 

shear-stress  gradient  and  is  assumed  to  he  1/2  A  ,  q  is  the  magnitude  of  the 

p 

velocity,  k  =  0.42  is  the  von  Karman  constant,  and  B  =  5.45.  Since  the  log- 
law  (IV- 32)  only  provides  the  velocity  magnitude,  in  order  to  obtain  the 
velocity  components,  it  is  assumed  that  the  velocity  vector  in  the  (x,y)  plane 
is  parallel  to  the  wall  and  in  the  (x,z)  plane  has  the  same  direction  as  at 
the  second  grid  point.  Since  the  second  grid  point  is  also  in  the  log-law 
region,  after  a  field  solution  has  been  obtained,  the  solution  at  the  second 
grid  point  can  be  used  to  update  the  guessed  value  of  and  the  procedure 
repeated  until  convergence.  Lastly,  conditions  on  p  and  p  are  not  required  on 
SB. 


Symmetry  Plane 


9_ 

9y 


A 

(U,W,p,p,k,e ) 


=  0 


(IV- 34) 


V  =  0 
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Symmetry  Planes  and  Sy 


(IV-35) 


W  =  0 


E.  Grid  Generation.  Grid  generation  is  an  important  aspect  of  the 
overall  numerical  procedures.  This  is  primarily  due  to  the  fact  that  the 
partially -parabolic  assumptions  are  made  with  reference  to  a  preselected  grid 
coordinate  curve  which  is  presumed  to  be  sufficiently  close  to  the  streamwise 
direction.  Also,  grids  that  are  not  sufficiently  smooth  can  cause  erroneous 
instabilities  in  the  solution.  Two  body-fitted  grid-generation  techniques 
have  been  used. 

For  the  small-domain  interactive  calculations  a  simple  analytic  grid- 
generation  technique  has  been  used  in  which  local  expanding  grids  for  each 
cross-section  are  pieced  together  in  the  streamwise  direction.  That  is,  the 
nonorthogonal  curvilinear  coordinates  (x,y,z)  are  defined  by 


x  =  X 


(IV-36) 


z  =  Az  y  Iq —  a  ^ 


z 
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where  Ay  -  y+v/U.  is  the  initial  y-direction  spacing,  y  is  the  expansion 

ratio  in  the  y-direction,  6(x,z)  is  the  houndary-layer  thickness,  Az  is  the 

initial  z-direction  spacing,  y  is  the  expansion  ratio  in  the  z-direction 

z  y 

and  a(x)  is  the  arclength  of  the  body  cross-section. 

For  the  large-domain  solutions  it  was  not  possible  to  use  such  a  simple 
analytic  technique  as  (IV— 36).  This  is  due  to  the  fact  that,  as  the  outer 
boundary  is  placed  at  further  distances  from  the  body  surface,  the  piecewise 
local-expansions  technique  produces  a  grid  in  which  the  coordinates  display 
humps  and  hollows  near  the  body  trailing  edge  and  midway  across  the  flow 
domain.  Such  irregularities  cause  instability  in  the  solution.  Consequently, 
for  the  large-domain  solutions,  a  more  sophisticated  body-fitted  grid-genera¬ 
tion  technique,  developed  by  Chen  and  Patel  (1985),  was  used  in  which  the 
coordinates  are  obtained  numerically  from  the  solution  of  a  set  of  Poisson 
equations  for  specified  boundaries  and  control  functions. 

F.  Global  Solution  Procedure.  In  the  previous  subsections  (IV.A)-(IV.E) 
the  various  aspects  of  the  viscous-flow-solution  method  have  been  described. 
The  inviscid-flow-solution  method  is  described  in  section  V.  In  this  section, 
the  global  and  interactive  solution  procedures  are  outlined.  The  main  steps 
are  as  follows: 

1)  Specify  all  boundary  and  initial  conditions,  including  an  initial 
guess  for  the  entire  pressure  field  and,  for  the  interaction  solu¬ 
tions,  the  boundary-layer  thickness. 

2)  Construct  the  grid  and  calculate  all  the  required  geometric  data. 
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3)  Evaluate  all  the  coefficients  in  the  momentum,  pressure-correction, 
and  pressure  equations. 

4)  Solve  the  momentum  equations.  An  under-relaxation  factor  ot  is  used. 

5)  Solve  the  pressure -correction  equation  and  correct  the  velocities. 

An  under-relaxation  factor  aA  is  used. 

P 

6)  Calculate  the  pseudovelocities  and  store  for  use  in  solving  the 
pressure  equation. 

7)  Evaluate  all  the  coefficients  in  the  k-e  equations  and  solve  the 

k-e  equations. 

8)  Repeat  steps  3)  -  7)  for  each  cross-section  until  the  exit  plane  is 
reached. 

9)  Calculate  the  displacement  body  and  for  the  interaction  solution  the 
new  inviscid-flow  solution. 

10)  Solve  the  pressure  equation  elliptically.  For  the  interaction  solu¬ 
tion  the  edge  conditions  are  updated  based  on  the  new  inviscid-flow 

solution.  An  under-relaxation  factor  a  is  used. 

P 

11)  Repeat  steps  3)  -  10)  until  convergence  is  reached  within  a  specified 
tolerance.  Actually,  the  convergence  criterion  used  for  the  results 
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to  be  presented  was  simply  that  there  was  no  noticeable  change  in  the 
body  surface  and  wake  centerline  pressure  distribution  when  viewed  on 
the  plotting  device.  Below,  the  notation  IT  is  used  to  refer  to  the 
global  iteration  number. 

V.  INVISCID  FLOW 

A.  Source-Panel  Method.  In  the  inviscid-flow  region  (region  1  in  figure 
1),  the  flow  is  assumed  to  be  irrotational.  A  velocity  potential  $  is  defined 
such  that  V$  is  the  perturbation  velocity  field,  i.e. 

A 

V  =  U  i  +  V$  (V-l) 

_E  0 

The  perturbation  potential  must  satisfy  the  Laplace  equation 

V2$  =  0  ( V-2 ) 


subject  to  the  boundary  condition 

$  =  -  D  n,  on  Sn  (V-3) 

n  o  l  d 

on  the  body, and  the  condition 

V4>  o  (V-4) 

at  infinity.  The  boundary-value  problem  (V-2)  -  (V-4)  for  the  perturbation- 
potential  can  be  solved  by  the  source -distribution  method.  The  conforming 
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panel,  source-panel  method  of  vou  Kerczek  et  al.  (1983)  has  been  used  here. 
The  perturbation  potential  is  expressed  by 


*  =  J  o  GdS 
SB 


(V-5) 


where  o  is  the  source  strength  and  G  is  the  Green  function 


with  X  =  (x,y,z)  the  field  point  and  =  (5  ,u  ,C  )  the  source  point.  Note  that 
equation  (V-5)  automatically  satisfies  condition  (V-4).  Substitution  of  (V-5) 
into  (V-3)  results  in  the  integral  equation  for  the  source  strength  cr,  i.e. 


( V— 7 ) 


Equation  (V-7)  is  solved  by  discretizing  the  body  surface  into  a  number  of 
conforming  surface  panels,  on  each  of  which  the  source  strength  is  assumed 
constant,  and  evaluating  the  integral  in  (V-7)  over  each  panel  approximately 
by  using  Gauss  quadrature. 

This  procedure  is  facilitated  through  the  use  of  a  parametric  surface 
equation  representation  for  the  body  surface: 


Rg  =  X  i  +  Y  (X,ct )  j  +  Z  (X,o )  k 


(V-8) 
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where  Rg  is  a  position  vector  which  describes  the  body  surface.  Here  X  is  the 
usual  lengthwise  coordinate  and  a  is  a  parameter  which  varies  along  the  girth 
of  each  section.  The  following  procedures  are  implemented.  First,  the  body 
offsets  for  one  side  are  mapped  into  the  X-<*  plane.  That  is,  corresponding 
(X,a)  pairs  are  assigned  for  each  given  offset  Rg.  Second,  the  derivatives, 
R^,  R^,  and  R^ ,  are  obtained  by  three-point  finite-difference  approxima¬ 
tions.  Third,  the  offsets  and  derivatives  are  interpolated  using  cubic- 
Hermite  splines  onto  a  uniform  {X^,  a^}  grid  chosen  so  as  to  provide  a  good 
surface  coverage  with  regard  to  surface  features.  With  the  surface  vector, 
Rg,  and  its  derivatives,  Rx,  R^  and  R^ ,  known  at  each  gridpoint,  their  values 
at  any  arbitrary  point  (X,a)  on  the  grid  can  easily  be  obtained  by  reinterpo¬ 
lation,  again  using  cubic-Hermite  splines.  Lastly,  (X,a )  are  specified  to 
form  the  input  needed  for  the  evaluation  of  the  integral  in  (V-7).  That  is, 
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(V-9) 


where  wk  is  the  weighting  factor  and  (Xk,  «k)  are  the  nodes  of  the  Gauss 
quadrature  formula.  The  resulting  system  of  linear  equations  is  solved  for  a 
by  Gauss-Seidel  iteration.  With  a  known,  the  surface  and  field  point  veloci¬ 
ties  are  readily  calculated,  again  using  Gauss  quadrature,  from 


V  =  U  i  +  /  a  VG  dS 
— a.  °  q 

and  the  pressure  is  obtained  from  the  Bernoulli  equation 
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B.  Displacement  Body.  As  explained  previously,  in  the  interaction 
calculation,  the  body  boundary  condition  (V-3)  is  applied  not  on  the  actual- 
body  surface  but  on  the  displacement-body  surface.  In  terms  of  the  present 
nonorthogonal  curvilinear  coordinate  system  the  displacement  body  is  defined 
by 


(V-10) 


where  V^  is  the  velocity  vector  of  the  outer  potential  flow  analytically 
continued  into  the  viscous-flow  region,  _V  is  the  viscous-flow  velocity  vector, 
#  and  A^  are  the  cross-sectional  areas  between  the  actual-body  surface  and 
the  displacement-body  surface  and  the  boundary-layer  surface  respectively, 

A  A 

and  dA  =  x  l^h^dydz  .  In  evaluating  (V-10),  for  axisymmetric  flow 

applications,  the  approximation  was  used  that  V  is  constant  across  A  *  i.e. 

P  5 


* 

6 


-  V)  •  i  b^h^  dy 


I  ^2h3  ^ 


(V-ll) 


This  approximation  can  be  removed  by  analytically  continuing  the  displacement- 
body  potential-flow  solution  inside  the  displacement  body;  however,  this 
requires  an  inviscid-flow  solution  method  based  on  the  symmetric  form  of 
Green's  theorem  (i.e.  source  and  dipole  distributions)  which  is  continuous 
across  the  body  surface  whereas  the  present  source-panel  method  is  not. 
Lastly,  some  comments  should  be  made  with  regard  to  the  evaluation  of  (V-10) 
for  three-dimensional  flow  applications.  A  straightforward  but  approximate 


34 


and  make  the  addition- 


procedure  is  to  again  assume  V-  is  constant  across  A, 
al  approximation  that  5*(x,z)  can  be  defined  in  terms  of  a  local  flux  balance, 

i.e. 


6* 

/  h^dy 


1 


(S  *  )•  idz 


/  (V  -  V)«  ih2h^dydz 


Landweber  (1986)  has  pointed  out  that  such  an  approximation  does  not  guarantee 
that  <5  *  is  a  stream  surface  of  the  continued  potential  flow  and  suggests  an 
alternative  method  for  evaluating  (V-10)  subject  to  an  explicit  condition 
that  6*  is  a  stream  surface. 


C.  Equivalent-Source  Method.  For  the  flat-plate  boundary-layer  and  wake 
test  case  (see  Section  VI. A.)  the  displacement  effect  of  the  boundary  layer 
was  included  through  the  use  of  the  equivalent -source  method.  In  this  case, 
such  a  method  is  acceptable  since  the  body  surface  is  flat.  A  two-dimensional 
potential  function  is  defined  by 


where 


X, 


$  =  U  X  +  /  aGdX 
o  ; 

*1 


G  =  £nR 


(V-12) 
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VI.  APPLICATIONS  FOR  TWO-DIMENSIONAL  AND  AXISYMMETRIC  FLOWS 


Results  are  now  presented  from  the  application  of  the  computational 
procedures  described  in  Sections  II  -  V  for  two-dimensional  and  axisymmetric 
flows.  The  two-dimensional-flow  application  is  a  simple  flat-plate  boundary 
layer  and  wake.  Results  are  presented  for  both  laminar  and  turbulent  flow. 
For  the  axisymmetric-flow  applications,  results  are  presented  for  turbulent 
flow  only,  and  for  two  of  the  family  of  afterbodies  for  which  experimental 
data  have  been  obtained  by  Huang  and  associates  at  DTNSRDC .  The  two  after¬ 
bodies  investigated  represent  examples  of  medium  and  strong  viscous-inviscid 
interaction.  In  the  discussions  to  follow,  all  coordinates  are  nondimension— 
alized  using  the  body  length  L,  with  X  =  0  at  the  body  leading  edge  (nose), 
and  velocities  and  pressure  are  normalized  using  the  free-stream  velocity  UQ 
and  the  fluid  density. 

A.  Flat  Plate  Boundary  Layer  and  Wake.  The  simple  case  of  a  flat-plate 

boundary  layer  and  wake  has  been  the  subject  of  many  previous  investigations. 
For  laminar  flow,  solutions  have  been  obtained  using  a  variety  of  approaches: 
thin-boundary-layer;  thin-boundary-layer  including  viscous-inviscid  interac¬ 
tion;  triple-deck  theory;  partially -parabolic  Navier-Stokes ;  and  Navier- 
Stokes.  Recently,  Chen  and  Patel  (1986)  have  performed  a  comprehensive  inves¬ 
tigation  in  order  to  establish  the  capabilities  of  their  partially -parabolic 
method,  extended  for  complete  Navier-Stokes  solutions,  by  comparing  their 
results  with  the  solutions  obtained  using  alternative  approaches  and  systema¬ 
tically  studying  the  influence  of  the  boundary  conditions,  the  size  of  the 
solution  domain,  and  the  grid  resolution.  A  direct  comparison  will  be  made 
between  results  using  the  present  approach  and  that  of  Chen  and  Patel.  For 
turbulent  flow,  results  have  also  been  obtained  using  a  variety  of  approaches; 
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however,  in  this  case,  a  direct  conparison  is  made  more  difficult  by  the  fact 
that  the  results  additionally  depend  upon  the  turbulence  model  utilized  and 
the  numerical  details  of  its  implementation.  Patel  and  Chen  (1986)  have  also 
performed  turbulent  flow  calculations  using  the  complete  Reynolds  equations 
and  the  k-e  turbulence  model  with  two  different  treatments  of  the  flow  close 
to  the  wall:  wall-function  approach;  and  an  eddy-viscosity  distribution.  The 
present  turbulent -flow  results  are  compared  with  both  the  results  of  Patel  and 
Chen  and  the  near-wake  experimental  data  of  Pot  (1979). 

i - laminar  Flow.  The  laminar-flow  calculations  were  performed  for  a 

UcoL 

plate  Reynolds  number  Rn  =  =  105  (where  L  is  the  plate  length)  which  is 
the  value  used  by  Chen  and  Patel  and  others.  Epical  large-  and  small-domain 
grids  used  for  the  calculations  are  shown  in  figures  5  and  6,  respectively. 
Referring  to  figures  5  and  1  for  notation:  the  number  of  axial  grid  points  is 
40;  the  number  of  transverse  grid  points  is  15;  the  inlet  boundary  Sj  is  at 
xu=.4;  the  exit  boundary  Sg  is  at  xd=2.5;  and  the  outer  boundary  SQ  is  at 
y0=12.  Referring  to  figures  6  and  1  for  notation:  the  number  of  axial  grid 
points  is  40;  the  number  of  transverse  grid  points  is  11;  the  inlet  and  exit 
boundaries  have  the  same  values  as  the  large-domain  grid;  and  the  outer  bound¬ 
ary  is  at  y0=1.26  6  where  6  (x)  is  the  boundary-layer  thickness  which  was 
specified  based  on  the  Blasius  solution.  The  y-direction  grid  expansion  for 
the  grids  shown  was  specified  based  on  the  turbulent  flow  conditions  since 
these  grids  were  actually  used  for  the  turbulent-flow  calculations  to  be 
discussed  subsequently.  The  y-direction  expansion  for  the  laminar-flow  calcu¬ 
lations  was  more  gradual.  The  Blasius  solution  was  used  to  specify  the  ini¬ 
tial  streamwise  profile  and  a  zero-gradient  condition  was  used  for  the  normal 
velocity  component  on  the  inlet  boundary  in  both  the  large-domain  and  small- 
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domain  calculations.  All  other  boundary  conditions  are  prescribed  as  dis¬ 
cussed  in  Section  IV. D. 

Figure  7  shows  the  pressure  distribution  on  the  surface  of  the  plate  and 
along  the  wake  centerline.  Results  are  shown  from  the  present  methods  both 
for  a  large  solution  domain  and  a  small  solution  domain,  including  viscous- 
inviscid  interaction  (interaction  solution).  Also  shown  for  comparison  are 
results  from  the  Navier-Stokes  calculations  of  Chen  and  Patel  (1986),  triple- 
deck  theory  (Melnik  and  Chow,  1975),  and  interactive  thin-boundary-layer 
theory  (Veldman,  1979).  It  is  seen  that  the  agreement  between  both  the  pre¬ 
sent  solutions  and  the  solution  of  Chen  and  Patel  for  x  =  .5421  is  excel- 

u 

lent.  The  reason  that  the  partially  parabolic  and  Navier-Stokes  solutions  are 
in  such  a  close  agreement  is  that,  at  this  high  Rn,  the  influence  of  stream- 
wise  diffusion  is  important  only  in  a  region  very  close  to  the  plate  trailing 
edge  and  not  resolvable  with  the  present  grid.  In  fact,  the  interactive  thin¬ 
boundary-layer  results  are  also  in  good  agreement,  indicating  that,  for  this 
very  simple  trailing-edge  flow,  the  influence  of  the  pressure  variation  within 
the  boundary  layer  is  also  small.  As  explained  by  Chen  and  Patel,  the  solu¬ 
tion  is  very  dependent  on  the  location  of  the  upstream  boundary.  Thus,  their 
solution  for  x  =  .1349,  which  has  also  been  included  on  figure  7  for  compari¬ 
son,  shows  large  differences.  Note  that  the  triple-deck  solution  for  which 

x  =  .1  is  consistent  with  the  solution  of  Chen  and  Patel  for  x  =  .1349  on 
u  u 

the  plate  but  not  in  the  wake.  In  the  wake,  the  triple-deck  solution  is  more 
consistent  with  the  other  solutions  obtained  using  larger  values  of  x^~  .5. 
Apparently,  this  is  due  to  the  treatment  of  the  downstream  boundary  in  the 
triple-deck  solution  which  is  matched  to  the  1/3  power  law  solution  at  too 
great  a  distance. 
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The  skin-friction  coefficient  and  the  wake-centerline  velocity  are  shown 
in  figures  8  and  9,  respectively.  As  seen  from  figure  8,  all  the  solutions 
for  the  skin-friction  coefficient  are  in  close  agreement  and  show  about  the 
same  deviation  from  the  Blasius  solution.  This  indicates  that  the  prediction 
of  the  skin-friction  coefficient  is  not  very  sensitive  to  the  details  of  the 
numerical  procedures  employed.  As  seen  from  figure  9,  the  differences  between 
the  solutions  for  the  wake-centerline  velocity  are  also  not  large  except  for 
the  triple-deck  solution  which  shows  larger  values  than  the  the  other  solu¬ 
tions  for  X  >  1.3  for  the  reason  discussed  previously.  The  close  agreement 
among  all  but  the  triple-deck  solutions  for  the  wake -centerline  velocity 
indicates  that  the  prediction  of  the  wake  properties  is  not  sensitive  to 
upstream  conditions  beyond  x^~  .5. 

Figures  10  and  11  show  the  converged  values  for  the  displacement  thick¬ 
ness  and  edge  velocity,  respectively,  used  in  the  interaction  solution. 
Referring  to  figures  10  and  11,  it  is  seen  that  the  magnitude  of  the  viscous- 
inviscid  interaction  for  the  flat-plate  boundary  layer  and  wake  is  weak; 
however,  its  features  are  characteristic  of  more  complex  trailing-edge 
flows.  Note  that  for  the  flat  plate  test  case  the  equivalent -source  method  is 
used  (see  Section  V. C.)  to  represent  the  displacement  effect  of  the  boundary 
layer  in  the  interaction  solution. 

The  interaction  calculations  were  started  with  free-stream  edge  condi¬ 
tions  (U  =  1. ,  and  pe=  0. )  which  were  subsequently  updated  after  each  global 
iteration.  The  large-domain  solution  converged  in  60  global  iterations  and 
the  interaction  solution  in  40.  In  both  solutions,  the  relaxation  factors 
used  were  (a,  a  ,  a*)  =  (.5,  .7,  1.).  Also,  both  solutions  required  about  5 

Jr  Jr 

minutes  cpu  on  a  Prime  9950  minicomputer. 
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2.  Turbulent  Flow.  The  turbulent  flow  calculations  were  performed  for  a 
plate  Reynolds  number  Rn  =  1.2  x  lO^1  which  results  in  a  momentum- thickness 
Reynolds  number  at  the  plate  trailing  edge  of  Rn0  =  3000.  This  corresponds 
fairly  closely  to  the  experimental  condition  of  the  wake  measurements  of  Pot 
(1979)  for  which  Rn0  =  2940.  The  large-  and  small-domain  grids  used  in  the 
calculations  are  shown  in  figures  5  and  6  respectively.  Referring  to  figures 
5  and  1  for  notation:  the  number  of  axial  grid  points  is  56;  the  number  of 
transverse  grid  points  is  15;  the  inlet  boundary  Sj  is  at  xu  =  .4;  the  exit 
boundary  SE  is  at  xd  =  16.25;  the  outer  boundary  SQ  is  at  yQ  =  1;  and  the 
first  grid  point  off  the  body  surface  was  located  at  y-*~40  .  Referring  to 
figures  6  and  1  for  notation:  the  number  of  axial  grid  points  is  56;  the 

number  of  transverse  grid  points  is  15;  the  inlet  and  exit  boundaries  as  well 
as  the  first  grid  point  off  the  body  surface  have  the  same  values  as  the 

large-domain  grid;  and  the  outer  boundary  is  at  yQ  =  1.26  6  for  X  <  1.  and 
varied  linearly  to  yQ  =  .25  at  the  exit. 

Figures  12  and  13  show  the  pressure  distribution  on  the  surface  of  the 
plate  and  along  the  wake  centerline  and  the  skin-friction  coefficient  respec¬ 
tively.  It  is  seen  that  the  large-domain  and  interaction  solutions  are  in 

excellent  agreement.  By  comparing  figures  7  and  12  it  is  seen  that,  for 

turbulent  flow,  the  extent  of  the  region  of  pressure  variation  is  reduced  for 
the  boundary-layer  region  upstream  of  the  trailing  edge  and  increased  for  the 
wake  region  downstream  of  the  trailing  edge.  Also,  the  pressure  recovery  in 
the  wake  occurs  with  a  favorable  pressure  gradient  for  laminar  flow  and  an 
adverse  pressure  gradient  for  turbulent  flow.  Also  shown  for  comparison  on 
figures  12  and  13  are  the  results  from  Patel  and  Chen  (1986)  using  the  k-e 
turbulence  model  with  wall  functions.  Note  that  their  results  are  for  a 
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slightly  different  Reynolds  number  (Rn  =  2.48  x  106)  and  that  their  grid 
resolution  was  similar  to  the  present  one  except  the  first  grid  point  off  the 
body  surface  was  located  at  y+~100  .  The  comparison  between  the  pressure 
distributions  is  excellent  for  X  >.9.  For  X  <.9  the  solution  of  Patel  and 
Chen  shows  lower  pressures.  This  is  due  to  the  initial  conditions  used  in 
their  calculations.  In  order  to  compare  skin-friction  coefficients  the  re¬ 
sults  of  Patel  and  Chen  were  Rn  scaled.  Referring  to  figure  13,  it  is  seen 
that  although  the  trends  are  identical  their  result  is  slightly  lower  than  the 
present  one.  This  is  due  to  the  larger  value  of  y+  for  the  first  grid  point 
off  the  body  surface  used  by  Patel  and  Chen. 

Figure  14  shows  the  usual  overall  parameters  for  describing  the  near 
wake,  the  half -width  b,  the  centerline  velocity  defect  wQ  =  1-UCL,  and  also 
the  shape  parameter  H  =  <5*/0.  Also  shown  for  comparison  is  Pot's  experimental 
data  and  the  results  of  Patel  and  Chen  for  wQ.  It  is  seen  that  the  agreement 
between  the  large-domain  and  interaction  solutions  for  w0  and  H  is  very  good; 
however,  both  solutions  deviate  from  the  experimental  data  for  100  _<  X/0  < 
600.  A  similar  deviation  from  the  experimental  data  is  seen  for  the  large- 
domain  solution  for  b.  The  interaction  solution  for  b  agrees  with  the  large- 
domain  solution  only  In  the  very  near  wake  and  elsewhere  shows  larger  val¬ 
ues.  In  general,  the  calculations  show  larger  values  for  b,  wQ  and  H  indicat¬ 
ing  a  thicker  wake  region.  The  agreement  between  the  present  results  for  w 
and  Patel  and  Chen  is  quite  good  with  the  differences  being  attributable  to 
the  different  values  of  y+  used  for  the  first  grid  point  off  the  body  sur¬ 
face.  Patel  and  Scheuerer  (1982)  have  compared  results  from  wake  calculations 
usipg  thin-boundary-layer  equations  and  the  k-e  turbulence  model  with  these 
same  experimental  data.  Their  results  show  better  agreement  in  the  near 
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wake  (X/0  <  350)  but  poorer  agreement  in  the  far  wake.  The  differences  be¬ 
tween  the  partially-parabolic  calculations  and  the  results  from  thin-boundary- 
layer  theory  are  consistent;  that  is,  the  influence  of  the  adverse  pressure 
gradient  in  the  near  wake  results  in  a  thicker  wake  region  in  the  partially- 

parabolic  solutions.  The  reason  for  the  poorer  agreement  of  the  partially- 

parabolic  solution  than  the  thin-boundary-layer  solution  with  the  experimental 
data  in  the  near  wake  is  not  known.  The  reason  for  the  discrepancy  between  the 
large-domain  and  interaction  solutions  for  b  is  due  to  the  width  of  the  inter¬ 
action  solution  domain  in  the  wake  (see  figure  6).  Further  calculations  using 
a  larger  growth  rate  for  5  in  the  wake  region  indicate  a  closer  agreement  of  b 
with  the  large-domain  solution. 

Figures  15  and  16  show  the  wake  momentum  thickness  and  centerline  eddy 
viscosity  respectively.  Referring  to  figure  15,  it  is  seen  that  the  present 
results  from  both  the  large-domain  and  interaction  solutions  indicate  larger 
values  of  0  than  the  experimental  data.  This  is  probably  due  to  the  differ¬ 
ence  in  the  trailing  edge  values  of  Rn0  .  Referring  to  figure  16,  it  is  seen 

that  the  present  results  indicate  a  larger  value  for  v  than  the  calculations 

"0 

of  Patel  and  Chen.  From  the  results  for  v  and  0  it  was  determined  that  the 

t 

asymptotic  value  of  v  /U  0  for  the  large-domain  solution  is  .0333  and  for  the 

t  o 

interaction  solution  is  .0406.  The  value  from  the  experimental  data  for  the 
range  400  <_  X/0  <_  1000  is  .035.  Patel  and  Scheuerer  and  Patel  and  Chen  ob¬ 
tained  values  of  .024  and  .022  respectively.  Figures  17  and  18  show  the 
asymptotic  (velocity-defect  w  =  1-U  and  stress  x  =  -uv  )  profiles  for  the 
large-domain  and  interaction  solutions.  Also  shown  for  comparison  is  Pot's 
experimental  data.  Both  solutions  for  the  velocity-defect  profile  show  excel¬ 
lent  agreement  with  the  experimental  data.  Both  solutions  for  the  stress 
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profile  are  in  agreement  and  show  a  peak  value  of  about  .029.  The  experimen¬ 
tal  stress  profile  indicates  larger  values  with  a  peak  value  of  about  .051. 
Patel  and  Scheuerer  and  Patel  and  Chen  also  obtained  lower  values  for  the 
stress  profile  than  the  experimental  data.  Their  peak  values  are  about  .033 
and  .03  respectively.  One  of  the  conclusions  of  Patel  and  Scheuerer  is  that 
the  k-e  model  does  not  adequately  predicb  the  observed  asymptotic  behavior. 
This  was  confirmed  in  the  present  investigation. 

Figures  19  and  20  show  the  converged  values  for  the  displacement  thick¬ 
ness  and  edge  velocity,  respectively,  predicted  in  the  interaction  solution. 
By  comparing  these  figures  with  figures  10  and  11  it  is  seen  that  the  influ¬ 
ence  of  turbulence  is  to  reduce  the  extent  of  the  viscous-inviscid  interac¬ 
tion. 

The  interaction  calculations  were  started  with  edge  conditions  (Ue  =  1, 
pe  =  0. )  which  were  subsequently  updated  after  each  global  iteration.  The 
large-domain  solution  converged  in  40  global  iterations  and  the  interaction 
solution  in  25.  In  both  solutions,  the  relaxation  factors  used  were 

(a,  a  ,  a" )  =  (.5,  .45,  1.).  Also,  both  solutions  required  about  5  minutes  of 
P  P 

cpu  on  a  Prime  9950  minicomputer. 

B.  Axisymmetric  Bodies.  For  axisymmetric  flow,  calculations  have  been 
performed  for  DTNSRDC  afterbodies  1  (Huang  et  al. ,  1978)  and  5  (Huang  et  al. , 
1980).  Afterbody  1  is  the  parent  form  of  a  family  of  three-dimensional  bodies 
with  elliptical  cross-sections  (Huang  et  al. ,  1983).  Referring  to  figure  21, 
which  shows  a  comparison  of  afterbodies  1  and  5,  it  is  seen  that  both  bodies 
have  similar  length/diameter  ratios  but  very  different  tail  forms.  In  parti¬ 
cular,  afterbody  5  (Ls/D  =  2.04)  is  blunter  than  afterbody  1  (Ls/D  =  4.3)  and 
has  greater  curvature  with  an  inflection  point.  Note  that  afterbody  5  is 
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almost  as  blunt  as  afterbody  3  (Ls/D  =  1.5)  which  exhibited  a  small  region  of 
flow  separation  near  the  tail  (.92  X  <_  .97)  in  the  experiments.  As  will  be 
discussed  next,  these  differences  in  body  geometry  result  in  both  a  larger 
exbent  and  magnitude  of  viscous-inviscid  interaction  for  afterbody  5  than 
afterbody  1. 

The  calculations  were  performed  for  afterbody  1  for  a  body-length  Rey¬ 
nolds  number  Rn  =  6.6  x  10^  which  corresponds  to  the  experimental  condition. 
The  large-  and  small-domain  grids  used  in  the  calculations  are  shown  in  figure 
22  and  23  respectively.  Referring  to  figures  22  and  1  for  notation:  the 
number  of  axial  grid  points  is  60;  the  number  of  transverse  grid  points  is  19; 
the  inlet  boundary  Sj  is  at  x^  =  .5;  the  exit  boundary  Sg  is  at  x^  =  16.25; 
the  outer  boundary  is  at  yQ  =  .8137;  and  the  first  grid  point  off  the  body 
surface  was  located  in  the  range  100  <  y+  <  160  .  Referring  to  figure  23  and 
1  for  notation;  the  number  of  axial  grid  points  is  60;  the  number  of  trans¬ 
verse  grid  points  is  11;  the  inlet  and  exit  boundaries  have  the  same  values  as 
the  large-domain  grid;  the  outer  boundary  is  at  yQ  =6  where  6(x)  is  the 
boundary-layer  thickness  which  was  specified  based  on  the  experimental  re¬ 
sults;  and  the  first  grid  point  off  the  surface  was  located  in  the  range 
80  <  y+  <  130  .  The  y-direction  grid  expansion  was  specified  such  that  the 
first  two  grid  nodes  are  within  the  log-law  region.  Simple  turbulent  flat- 
plate  profiles  were  used  to  specify  the  initial  conditions  in  both  the  large 
and  small  domain  solutions.  All  other  boundary  conditions  are  prescribed  as 
discussed  in  Section  IV. D. 

Figure  24  shows  the  pressure  distribution  on  the  surface  of  the  body  and 
along  the  wake  centerline.  Results  are  shown  from  the  present  methods  both 
for  a  large  solution  domain  and  a  small  solution  domain,  including  viscous-in- 
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viscid  interaction.  Also  shown  for  comparison  are  the  experimental  results 
and  the  inviscid-flow  solution  without  interaction.  It  is  seen  that  both  of 
the  present  methods  are  in  good  agreement  with  the  experimental  results. 
Actually,  the  interaction  solution  appears  to  be  slightly  in  better  agreement. 
This  is  probably  due  to  the  influence  of  the  initial  conditions  which,  in  the 
case  of  the  interaction  solution,  include  a  more  proper  matching  with  the 
external  flow.  Also,  the  grid  resolution  within  the  boundary-layer  region  is 
better  for  the  interaction  solution.  Note  the  large  gradients  in  pressure 
exhibited  in  both  solutions  in  the  immediate  vicinity  of  the  trailing  edge  (X 
=  1).  A  part  of  this  behavior  is  no  doubt  a  result  of  the  rapid  change  in 
curvature  of  the  streamlines  associated  with  the  closing  of  the  body  and 
transition  into  the  wake.  However,  it  was  also  found  that  the  solution  in 
this  vicinity  is  sensitive  to  the  grid  and  detailed  numerical  treatments  at 
the  trailing  edge.  A  comparison  of  the  present  results  with  the  inviscid-flow 
solution  without  interaction  provides  one  indication  of  the  magnitude  of  the 
viscous-inviscid  interaction  (47#  reduction  in  the  maximum  value  of  Cp  at  x  = 
.975).  This  comparison  is  made  somewhat  difficult  by  the  inaccuracy  of  the 
inviscid-flow  solution  very  near  the  trailing  edge.  In  the  present  inviscid- 
flow  method,  as  is  the  case  with  most  other  singularity-distribution  methods, 
the  solutions  are  not  accurate  in  regions  where  the  angles  between  adjacent 
panels  are  small.  This  inaccuracy  is  well  known  and  was  not  removed  since  the 
inviscid-flow  solution  without  interaction  is  not  used  in  the  present  viscous- 
inviscid  interaction  approach. 

Figure  25  shows  the  wall-shear  velocity  U  and  is  of  similar  format  as 
the  previous  one.  The  experimental  values  were  obtained  from  Preston  tube 
measurements.  In  this  case,  the  comparison  between  the  calculations  and  the 
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experimental  values  is  not  quite  as  good.  For  .9  <  X  <  1.  the  calculations 
show  larger  values  than  the  experiment.  This  discrepancy  is  attributed  to  the 
well  known  deficiencies  of  the  k-e  turbulence  model  for  thick-boundary-layer 
flow  as  will  be  discussed  further  subsequently.  Also,  the  large-domain  solu¬ 
tion  predicts  slightly  higher  values  for  than  the  interaction  solution  for 
X  <  *9.  This  is  no  doubt  due  to  the  influence  of  the  initial  conditions. 

Figure  26  shows  the  wake-centerline  velocity  U  and  is  of  similar  format 

GLi 

as  the  previous  ones.  The  agreement  between  both  the  calculation  methods  and 
the  limited  experimental  data  is  again  quite  good;  however,  the  calculations 
show  a  slower  recovery  than  that  indicated  by  the  limited  near-wake  experimen¬ 
tal  data.  The  largest  differences  between  the  large  domain  and  interaction 
solutions  are  in  the  intermediate -wake  region  2.  <  X  <  5.  and  are  consistent 
with  differences  in  pressure  as  shown  in  figure  24. 

Figures  27  and  28  show  the  convergence  history  (pressure  distribution  on 
the  surface  of  the  body  and  along  the  wake  centerline  and  the  displacement 
thickness)  for  the  large-domain  and  interaction  solutions  respectively. 
Values  are  shown  for  every  five  global  iterations.  The  planar  definition  of 
displacement  thickness  has  been  used  for  the  large-domain  solution.  The 
interaction  calculations  were  started  with  free-stream  edge  conditions  (Ue  = 
1.,  pe  =  0.).  After  twenty  iterations,  the  edge  conditions  were  updated  using 
the  latest  value  of  displacement  thickness.  Subsequently,  the  edge  conditions 
were  updated  every  five  global  iterations  until  convergence  was  achieved.  A 
comparison  of  figures  27  and  28  shows  that  the  convergence  characteristics  of 
the  two  solutions  are  quite  different.  The  large-domain  solution  converges 
monotonically  in  50  iterations.  The  interaction  solutions  converge  with 
oscillations  in  40  iterations.  Basically  the  interaction  solution  converges 


46 


in  two  stages.  The  first  stage  is  with  free-stream  edge  conditions  and  leads 
to  an  underprediction  of  both  Cp  and  6  .  The  second  stage  is  with  the  dis¬ 
placement-body  edge  conditions  and  the  solution  converges  quite  rapidly. 

29  shows  the  iteration  history  of  the  edge  velocity.  The  small  changes 
in  the  displacement-body  shape  after  20  global  iterations  lead  to  impercep¬ 
tible  changes  in  Ue.  Figure  30  shows  a  comparison  between  the  converged  edge 
pressure  obtained  from  the  displacement  body  with  that  obtained  from  the 
actual  body,  and  the  actual  body  including  the  equivalent-source  method  to 
represent  the  displacement  effect  of  the  boundary  layer.  It  is  seen  that  the 
displacement-body  edge-pressure  maximum  is  shifted  upstream  as  compared  to  the 
actual-body  edge  pressure  which  results  in  greater  edge  pressures  for  X  <  .95 

/V 

and  lower  values  for  X>  .95.  Note  that  the  equivalent-source  method  edge 
pressure  shows  the  correct  tendency  but  with  only  a  minimal  modification  to 
the  actual-body  result.  This  clearly  demonstrates,  as  was  discussed  pre¬ 
viously,  that  the  equivalent-source  method  is  inaccurate  for  bodies  with 
noncusped  trailing  edges.  Figure  31  shows  a  comparison  of  the  experimental 
and  interaction  solution  displacement  thickness.  The  calculated  values  are 
slightly  below  the  experimental  values. 

Lastly,  for  afterbody  1,  figure  32  shows  the  solution  profiles 

(U>v>P>L,e)  for  a  number  of  X-stations  between  the  inlet  and  the  outlet. 

Y— R 

Note  that  in  these  figures  the  radial  coordinate  has  been  defined  as  Y  =  ^ _ 9 

max 

where  RQ(x)  is  the  local  body  radius  and  is  the  maximum  body  radius. 

Wherever  possible,  a  comparison  has  been  made  with  the  experimental  data.  In 
general  the  agreement  between  the  large-domain  and  interaction  solutions  is 
very  good  and  consistent  with  the  previous  discussions.  Also,  both  solutions 
show  good  agreement  with  the  experimental  data.  However,  note  the  fact  that 
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while  both  solutions  indicate  similar  transverse  pressure  gradients  Py  there 
is  a  systematic  difference  in  pressure  magnitude.  The  large-domain  solution 
predicts  lower  pressures  than  the  interaction  solution  which  in  general  shows 
better  agreement  with  the  experimental  data.  The  reason  for  this  difference 
is  not  known.  It  may  be  related  to  the  outer  boundary  conditions  in  the 
large-domain  solution  since  its  characteristics  are  similar  to  a  blockage 
effect.  The  principal  differences  between  the  calculations  and  the  experi¬ 
ments  is  that  the  calculations  tend  to  overpredict  the  velocity  and  turbulent 
kinetic-energy  profiles.  This  result  is  due  to  deficiencies  of  the  standard 
turbulence  model  which  is  known  to  overpredict  the  level  of  turbulence  in 
thick  boundary  layers,  presumably  caused  by  the  use  of  an  isotropic  eddy 
viscosity  and  the  neglection  of  curvature  effects.  Another  cause  may  be  the 
use  of  wall  functions.  Note  the  significant  variation  in  pressure  across  the 
boundary  layer  in  the  vicinity  of  the  trailing  edge  .95  <  X  _<  1.05  indicating 
the  necessity  of  including  such  effects  in  modeling  the  present  flow. 

The  calculations  were  performed  for  afterbody  5  for  a  body-length  Rey¬ 
nolds  number  Rn  =  9.3  x  10^,  which  corresponds  to  the  experimental  condi¬ 
tion.  The  large-  and  small-domain  grids  used  in  the  calculations  are  shown  in 
figures  33  and  34  respectively.  Referring  to  figures  33  and  1  for  notation: 
the  number  of  axial  grid  points  is  60;  the  number  of  transverse  grid  points  is 
19;  the  inlet  boundary  Sj  is  at  x  =  .5;  the  exit  boundary  Sg  is  at  x^  =  16.25; 
the  outer  boundary  is  at  yQ  =  .8137;  and  the  first  grid  point  off  the  body 
surface  was  located  in  the  range  120  <  y+  <  200  .  Referring  to  figures  34  and 
1  for  notation:  the  number  of  axial  grid  points  is  60;  the  number  of  trans¬ 
verse  grid  points  is  11;  the  inlet  and  exit  boundaries  have  the  same  values  as 
the  large-domain  grid;  the  outer  boundary  is  at  yQ  =  5  where  5  (x)  is  the 
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boundary  layer  thickness  which  was  specified  based  on  the  experimental  re¬ 
sults;  and  the  first  grid  point  off  the  body  surface  was  located  in  the  range 
100  <  y+  <  160  .  The  y-direction  grid  expansion  was  specified  such  that  the 
first  two  grid  points  are  within  the  log-law  region.  Simple  turbulent  flat- 
plate  profiles  were  used  again  to  specify  the  initial  conditions  in  both  the 
large-  and  small-domain  solutions.  All  other  boundary  conditions  are  pre¬ 
scribed  as  discussed  in  Section  IV.  D.  In  the  discussion  to  follow,  all 
figures  for  afterbody  5  are  of  a  similar  format  as  that  described  earlier  for 
afterbody  1. 

Figure  35  shows  the  pressure  distribution  on  the  surface  of  the  body  and 
along  the  wake  centerline.  Comparing  figures  35  and  24,  it  is  seen  that  the 
pressure  distribution  in  the  tail  region  of  afterbody  5  shows  even  more  radi¬ 
cal  variations  than  on  afterbody  1.  This  includes  both  a  lower  minimum  pres¬ 
sure  upstream  of  the  trailing  edge  and  a  higher  maximum  pressure  at  the  trail¬ 
ing  edge.  It  is  seen  that  both  the  large-domain  and  interaction  solutions  are 
in  good  agreement  with  each  other  and  the  experimental  data.  The  level  of 
agreement  is  about  the  same  as  that  obtained  for  afterbody  1.  The  magnitude 
of  the  viscous-inviscid  interaction  is  larger  for  afterbody  5  than  for  after¬ 
body  1  (50%  reduction  in  the  maximum  value  of  Cp  at  X  =  .975). 

Figure  36  shows  the  wall  shear  velocity  and  figure  37  the  wake  center- 
line  velocity  Dct  The  level  of  agreement  is  not  as  good  as  that  obtained  for 
afterbody  1.  Figure  38  shows  the  converged  value  of  displacement  thickness 
and  figure  39  the  resulting  edge  velocity.  A  comparison  of  figures  38  and  28 
shows  that  the  viscous-inviscid  interaction  is  larger  for  afterbody  5  than 
afterbody  1,  resulting  in  an  increased  displacement  thickness  for  afterbody 
5.  The  displacement  thickness  is  also  shown  in  figure  40  where  it  is  conpared 
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with  the  experimental  data.  The  agreement  is  very  good.  Lastly,  the  detailed 
solution  profiles  at  various  X-stations  are  shown  in  figure  41.  The  results 
are  similar  and  consistent  with  those  described  earlier  and  do  not  require 
further  elaboration.  However,  note  the  increase  in  pressure  variation  across 
the  boundary  layer  as  compared  to  afterbody  1,  again  indicating  the  larger 
viscous -inviscid  interaction  and  thick-boundary-layer  effects  for  afterbody  5 
than  1. 

For  both  afterbodies  1  and  5,  the  large-domain  solutions  converged  in  50 
global  iterations  and  the  interaction  solutions  in  40.  For  afterbody  1,  the 
relaxation  factors  were  (a,  a  cT)  =  (1.,  .2-. 5,  1.)  for  the  large-domain 

-H  Jr 

solution  and  (a,  a^,  ct^)  =  (.6,  .2-. 5,  1.)  for  the  interaction  solution.  For 

afterbody  5,  the  relaxation  factors  were  (a ,  a  ,  a*  )  =  (.6,  .2-. 5,1.)  for  the 

P  P 

large-domain  solution  and  (a,  ct^,  a*  )  =  (.5,  .2-. 5,  1.)  for  the  interaction 

solution.  In  all  cases,  in  the  wake  a  =  .05  for  IT  <  10  and  a  =  .1  -  .2  for 

P  P 

IT  >_  10.  Also,  in  all  cases,  the  solutions  required  about  10  minutes  of  cpu 
on  a  Prime  9950  miniconputer. 

VII.  CONCLUDING  REMARKS 

It  has  been  shown  that  trailing-edge  flows  with  thick  boundary  layers  can 
be  modeled  using  viscous-inviscid  interaction  procedures  if  all  the  important 
aspects  of  the  flow,  namely  the  pressure  variation  across  the  boundary  layer 
and  the  displacement  effect  of  the  viscous  flow  on  the  external  inviscid  flow, 
are  taken  into  account.  The  validation  of  the  present  interaction  procedures 
has  been  accomplished  by  comparing  the  results  for  two-dimensional  and  axisym- 
metric  flows  with  large-domain  solutions,  in  which  the  entire  zone  of  viscous- 
inviscid  interaction  is  captured,  obtained  using  the  same  numerical  proced- 
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ures.  Also,  a  direct  comparison  has  been  made  with  other  methods,  including 
the  finite -analytic  method  of  Chen  and  Patel  (1985),  and  with  available  exper¬ 
imental  data. 

In  general,  the  comparison  between  the  present  interaction  and  large- 
domain  solutions  is  very  satisfactory  for  all  the  cases  investigated.  Some 
small  differences  are  evident,  as  might  be  expected.  The  direct  comparison 
between  the  present  methods  and  other  methods  and  also  the  experimental  data 
shows  excellent  agreement.  Of  particular  interest  is  the  comparison  with  the 
method  of  Chen  and  Patel;  since,  their  method  and  the  present  ones  have  cer¬ 
tain  features  in  common  and  are  quite  different  in  other  respects.  Specifi¬ 
cally,  the  velocity-pressure  coupling  procedures  as  well  as  the  turbulence 
model  are  identical;  however,  the  coordinate  systems  used  in  solving  the 
governing  equations,  the  discretization  procedures  employed,  as  well  as  other 
numerical  treatments  are  very  different.  Chen  and  Patel  (1985)  and  Patel  and 
Chen  (1985)  also  present  results  for  afterbodies  1  and  5,  including  a  compari¬ 
son  with  the  same  experimental  data.  The  level  of  agreement  is  very  similar 
to  that  shown  with  the  present  methods.  Thus,  it  would  seem  that  the  most 
critical  aspect  of  computational  methods  for  thick-boundary-layer  trailing- 
edge  flows  is  the  velocity-pressure  coupling  rather  than  the  discretization 
procedure  (finite  analytic  vs.  finite  difference)  and  other  numerical  proced¬ 
ures  employed. 

As  to  the  relative  advantages  of  the  interaction  vs.  large-domain  solu¬ 
tions,  the  latter  does  not  require  any  approximations  with  regard  bo  the 
viscous-inviscid  interaction  as  does  the  former;  however,  the  interaction 
solution  was  shown  to  be  just  as  accurate  as  the  large-domain  solution,  indi¬ 
cating  that  the  present  viscous-inviscid  interaction  procedures  can  capture 
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this  important  influence  on  the  flow.  With  regard  to  computational  effici¬ 
ency,  when  the  equivalent -source  method  is  used  the  interaction  solution  is 
more  efficient;  however,  when  the  displacement-body  method  is  used,  which  is 
the  case  for  bodies  of  interest,  then  the  large-domain  solution  is  more  effi¬ 
cient.  This  is  because  of  the  additional  computational  effort  in  calculating 
the  inviscid  flow.  The  prescription  of  the  inviscid  flow  at  the  boundary- 
layer  edge  does  speed  up  the  convergence  rate  of  the  viscous-flow  solution  but 
not  enough  to  offset  the  increase  pointed  out  above.  Note  that  we  have  used  a 
three-dimensional  source-panel  method  for  the  inviscid-flow  solution.  For 
three-dimensional  applications,  it  is  expected  that  the  interactive  solution 
will  be  more  attractive  since  the  saving  in  computational  effort  in  the  vis¬ 
cous-flow  solution  will  have  a  more  substantial  effect  on  the  overall  computa¬ 
tional  efficiency.  Calculations  for  three-dimensional  body  geometries  are  now 
in  progress  and  will  be  reported  on  in  the  future.  Lastly,  it  should  be 
mentioned  that,  besides  the  academic  interest  in  the  nature  of  an  interaction 
solution,  it  also  has  a  practical  value.  For  example,  the  modification  of  the 
external  inviscid-flow  solution  due  to  the  thick  boundary  layer  and  wake  is 
obtained  as  part  of  the  solution,  and  for  applications  to  ship  boundary  lay¬ 
ers,  it  is  readily  extendable  to  nonzero-Froude-number  calculations. 
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Figure  1.  Definition  Sketch  of 
Flow-Field  Regions 


Figure  2.  Nonorthogonal  Curvilinear 
Coordinates 
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Figure  3.  Finite-Difference  Molecule 
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Figure  4.  Definition  Sketch  of  Boundary  Surfaces 
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Figure  5.  Large-Domain  Grid  for  Flat-Plate  Boundary  Layer  and  Wake 
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Figure  6.  Small-Domain  Grid  for  Flat-Plate  Boundary  Layer  and  Wake 
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Figure  7.  Pressure  Distribution  on  the  Surface  of  the 
Plate  and  along  the  Wake  Centerline 


(C,  -s/RtT  )/2 


r 


1.0 


0.9 


0.8 


0.7 


06 


0.5 


04 


\ 


\ 


X 


\ 


X 


\ 


LAMINAR  FLOW 
•INTERACTION 
•LARGE  DOMAIN 


- BLASIUS  SOLUTION 

•  CHEN  AND  PATEL  (xQ=0.5421) 

□  CHEN  AND  PATEL  (xu  =  0.1349) 

O  TRIPLE-DECK  SOLUTION 

A  THIN  BOUNDARY  LAYER 
INTERACTION 


X 


X 


X, 


^r 


0.3 


0.0 


0.2 


04 


0  6 


0  8 


Figure  8.  Skin-friction  Coefficient 
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Figure  9.  Wake-Centerline  Velocity 
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Figure  10.  Displacement  Thickness 
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Figure  11.  Edge  Velocity 
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Figure  12.  Pressure  Distribution  on  the  Surface  of  the 
Plate  and  Along  the  Wake  Centerline 


Figure  13.  Skin-friction  Coefficient 
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Figure  14.  Near-Wake  Parameters 
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Figure  15.  Wake  Momentum  Thickness 


Figure  16.  Wake~Centerline  Eddy  Viscosity 
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Figure  17-  Large-Domain  Solution  Asymptotic  Profiles 
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Figure  18.  Interaction  Solution  Asymptotic  Profiles 
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Figure  20.  Edge  Velocity 


69 


max 


AFTERBODY  1 

AFTERBODY  5 

L 

3.066  m 

2.91  m 

D 

0.2794  m 

0.28  m 

Ls/L 

0.39 

0.196 

Ls/D 

4.3 

2.037 

L/D 

10.97 

10.39 

Figure  21.  Geometric  Data  for  DTNSRDC  Afterbodies 
1  and  5 


70 


0.05 


0.04 

0.03 

0.02 

0.01 


-  INTERACTION 

- LARGE  DOMAIN 

O  EXPERIMENT 


O 


0.00  1 - - - - - ^ ^ ^ ^ - 1 

0.6  0.7  0.8  0.9  1.0 

X 


Figure  25-  Wall-Shear  Velocity 
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Figure  26.  Wake-Centerline  Velocity 
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Pressure  Distribution  on  the  Surface  of  the  Body  and  Along 
the  Wake  Centerline 
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Figure  27.  Convergence  History  of  Large-Domain  Solution 
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Figure  28.  Convergence  History  of  Interaction  Solution 
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Figure  29.  Edge  Velocity 
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Figure  30.  Edge  Pressure 
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Figure  31.  Boundary  Layer  and  Displacement  Thicknesses 
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Figure  32.  Solution  Profiles 
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Figure  32.  (Continued) 
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Figure  32.  (Continued) 
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Figure  35.  Pressure  Distribution  on  the  Surface  of  the  Body 
and  Along  the  Wake  Centerline 
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Figure  37-  Wake-Centerline  Velocity 
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Figure  38.  Displacement  Thickness 
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Figure  39 .  Edge  Velocity 
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Figure  40.  Boundary  Layer  and  Displacement  Thicknesses 
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Figure  41.  Solution  Profiles 
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Figure  41.  (Continued) 
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APPENDIX  I:  Equations  in  Nonorthogonal  Curvilinear  Coordinates 


In  Section  IV. A  the  procedures  for  obtaining  the  governing  equations  in 
nonorthogonal  curvilinear  coordinates  have  been  discussed.  Since  the  result¬ 
ing  equations  are  quite  lengthy  they  are  provided  in  this  Appendix.  First, 
the  continuity  and  momentum  equations  are  presented.  Subsequently,  the  equa¬ 
tions  of  turbulent -kinetic -energy  k  and  its  dissipation-rate  e ,  which  are  used 
to  model  the  Reynolds  stresses  in  the  momentum  equations,  are  presented.  The 
partially-parabolic  forms  of  the  Reynolds  and  turbulence-model  equations  are 
provided  in  Appendix  II. 

A.  Continuity  Equation 


where 


with 


1  9U  19V  1  9W  TT 

h^  3l  +  h2  9y  +  h^  87  +  Ua  +  ^  +  Wc 


=  0 


a  = 


b  = 


c  = 


1 

8h2+  1 

3h3  1 

9S 

hlh2 

3x  h1h( 

9  x  h  S 

9  x 

1 

!lh  t  _u_ 

9hl+  1 

9  S 

h2h3 

ay 

9y  v,  Q 

3y 

1 

’V,  1 

3V  1 

9  S 

h2h3 

9  z  h^h^ 

9  z  h^S  9  z 

S  =  s/h^h^cij 


( A— 1 ) 


(A-2) 


(A- 3) 


and  s  is  the  triple  product  defined  by  (IV-13). 

B.  Reynolds  Equations 

1.  x-Momentum  Equation 


U  9  U  V  9  U  W  9  U  19  u^  1  9  uv  1  9  uw 
h^9x'  +  h^9y  +  h^9z'  +  h^  97"  +  h^  Ty  +  h^  7z 


+  R^UV  +  (R^+  b)  uv 


94 


+  R0(V2+  V2)  +  R-  (W2+  w2)  +  R.U2  +  (R.  +  a)u2  +  RC(VW  +  vw) 
^  i  4  4  5 


*  R6uw  ♦  (R6  *  c)  uw  =  -  [sin2*  *T  ^|f  ^|fl 

,_2IT  „  3U  „  .  3U  „  3U  3 V  3V  3V  3W  3W 

+  v  [V  J  +  C1  3x  +  C2  3y  +  C3  3z  +  C4  3x  +  C5  3y  +  C6  Tz  +  C7  Tx  +  C8  3y 


3  W 

*C9S*  C10U  *  C11V  +  C12W1 


(A-4) 


where  (a,b,c)  and  S  are  given  by  (A-2)  and  (A-3),  respectively,  and 


a  =  cos  v  cos  y  -  cos  X 


0  =  cos  X  cos  v  -  cos  y 


(A-5) 


y  =  cos  y  cos  X  -  cos  v 

with  (X  ,y  ,v  )  defined  by  (IV-11).  The  Laplacian  operator  in  nonorthogonal 
curvilinear  coordinates  is  defined  by  (where  Q  is  a  scalar  variable) 


w2n  _  A  A_  r»  i Q  x  3Q  £Q,  ^  1  3  r  3Q  3Q  3Qi 

7  Q  "  s  3x  [all  3x  a12  3y  a13  3z]  s  3y  ta21  3x  +  a22  3y  +  a23  3z] 


13,  3Q  3  Q  3  Q . 

s  3 z  [a31  3x  +  a32  3y  +  a33  3z] 


(A-6) 


where 


2  2  2  2 
h2h^sin  X  h-^h-Y 

all  =  I  >  a12  =  I 


hlh2h36 


13  s 
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*22 


,2,2  .  2 

h^h^sin  y 


23 


2 

h-^h^ 


,2,  2  .  2 
h^hgSin  v 


*33 


( A-7) 


and  -  ajj-  for  i  f  j.  The  Rj_  (i  -  1,6)  and  (i  =  1,  12)  coefficients 
are: 


3,1  a 


1  a 


*1  =  “2  [h3T  9^  (h3  cos  y  1  "  017  9^  (h2COS  V  1  " 


9h-, 

cosv  i  la 


S  2  3 


2  3 


a  z  i^h  a  x  u  3 
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hlh2S 


a  x 


W 


9y 


R. 


3 


a  (h  cos  A )  ah_  .  2,  3h„ 

,  3 _  _ 2,  sm  A  ,  2  a 


2  2 
h2h3S 


9y 


3z  v/  3x  ' ^  U‘1 


(h-,  cos  v  )  ] 


R3  = 


R,  = 


sin^A  .3 


3h. 


h^h^S' 


o  3h 

2  [£f  (hicos  p )  -  air1  -  ~j—  r~^  -  “ 


h2h3s 


2  9y  3  z  2 


ah. 


ah. 


h-^S' 


hlh2S 


ah. 


R5  = 


h2h3S 


2  3y 


-  cos  A 


9ho  •  K 
2 ,  sin  A 

az  J  _ 

hlh2S 


a  h 


(A-8) 


(A-9) 


(h„  cos  A ) ]  (A- 10) 


3  X  V  **-|  JN 

2  l37  (h3  003  11 1  "  a^1  -  7JLr  [ay  -  a7  <h2  cosv)l  U-lD 


2  a 


2  tcos  A  -  gy  (h1cos  y  )  ]  + 

ah„  ah„ 


te  lhi cos  v)  -  ai  (h3  cos  x)|  -rrT  lcosX  9/-32 


JL.hjS 


•  2, 
sm  A 

h-jh-^S* 


(A-12) 


.  2  ah  ah 

n  _  sm  A  ,  1  3, 

-  - —  lx—  -  cos  y  1 


6  2  l8  z 

h^h^S 


3  x 


3  3hl  3h3 

— 2  [cos  y  a^  -  rr1 


h^JS 


R2^3^ 


ah 

■hf  tcos  M  a y 


If  (h2  cos  v  }  1  +  7JLT2  [tlT  {h2cos  X)  ~h  (hlcos  W  )  1 


hlh2S 


(A-13) 


2,  .  2 


2.  2 


C  p  _  J_  ri_  hih2h3sin  x  a_  hih2h3y  a  /hih2h3e , , 

1  sl^  3x  s  3  y  s  3z  s  ^ 


2  2.2 

X  3  g  h-^hph^cos  v  ^ 

—3 - 57  <hT>  * - — -  ' 

s  Is 


ah. 


aF  (h2cos  x  1  "  aF1 
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h  h„h„cos  y  .  ahp  h.h^  Lh  cos  y 

*  - 5 -  '37  <V°S  -  j/1  *  4^  '=os  1  W  {—s - ’ 

h  h  cos  v  h  h  h  h  cos  v  .  h0h„cos  y 

-  57  ( — i - 1 1  +  1  lcos  x  37  < - i - 1  -  37  < - 7 - > 1 


(A-14) 


2h^h^cos  y  g 


g  h^l^h^cos  X 

(l^cos  v )  -  gy-  (h^  cos  y  )  ]  +  - 2 - 


g  3h^  h-^h^  g  g 

^3"x  ^h3cos  v  )  ”  TzT^  +  — 2~  ^Tz  (hih2cos  X  ^  “  G0S  ^  Yz  (hih2GOS  v  ^  1 

s 
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-  — 2  th^  gY  (h2cos  v  )  +  h-^  g-^-  (h2cos  X  )  ]  +  - 2 - gy  (h^h^cos  y  ) 


hlh2h3Y  3 _  .  s _ .  _  1  3_  hlh2h3Y 

3  3x  hn  s  3x  s 
s  1 


(A-15 ) 


h|h|h3p 

u3  3  3x  'h. 


hlh2 


3  s  ^1^2  3 
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+  — —  [cos  v  g-y  (h^cos  y  )  -  g^r  (h^cos  X)]  +  - -  [gy  (h^cos  y 


3  n  n  hlh2  .  2  9hl  3  13  hlh2h3?  v 

3^  (h2cos  v  )  ]  +  —2-  [cos  v  gj-  -  g^  (h3cos  y  )  ]  -  -  ^  ( ) 

s 


h  h  h  cos  X  3h 

JLST -  [h  (h2cos^  -3 3T1 


(A-16) 


2  2  2 
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n  d,Sx  1  J  r  ^  d  .  2.  x  3 
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4 


.3 


3y  h2'  s 


3y 
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-^0^0 


h,h  h„h  h„h  cos  X  h.h0h0cos  X 
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3#J  +  “Hr^  fc  (h2cos  x }  -  33T] 
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hjh^sin  X  9h2  h1h2h^Y  i  3  s  l  9h2  hlh33  9h2 
C5  "  _  2  3 x  2  s  3  y  (h0'  2  3y  1  2  3z 

S  S  c  S 

h  h  cos  u  h.  .  h  h  cos  X  h  3h  h  h  cos  X 

- - i -  [7“  37  (h3cos  x  5  +  37  ( - I - )  -  “  37“] - I - 

h2  3  h2  3  hlh3  ,h3  9h2  h3  3  ,h  \ 

[“37  (hlCOS  v)  -  —  37  (h3COS  X)1  "  “7“  [“7T  "  “37  (hlCOS  V) 


9  h,h~COS  V  h,h,,  3h„  3 

-  ^7  ( - - - )  ]  +  — 2  [cos  v  37“  ~  cos  x  (hQcos  \i  )  ] 


3y 


hlh2 

+  —  loos  v  ^  (  s 


3  ,hlh2°0S  X 


,  ,  3  7*1  3 

)  -  cos  X  —  ( - “ 

9  Z  S 


3  y  Vil3" 
h,  h„cos  v 


-)] 


hlh3  3  3 

+  — —  [cos  X  (h2cos  v)  -  cos  v  g-^  (h2cos  X)] 

s 


(A-l  8) 


2  2 

h,  h0h03  „  h,h„  .  hnh„  .  h.,hncos  v 

°6  = - 5 - 3y  (h7  ~s  [oos  11  37  (~>  -  008  X  57  ( - 5 - 

S  <L 


)] 


2  2 
h,h2  „  a  h,h2cos  v  3h2  „ 

+  — 2~  [gy  (h^cos  U )  -  (h2cos  V  )  ]  +  - c— -  [g-^-  -  gj  (h^cos  X  )  ] 

s  s 

L  hlh2COS  V  rhl  3  hl  3h2  3  ,hlh2,,  hlh2  fh2  3  „ 

+  - I -  [“  37  (h3C°S  X  }  “  “  3^  “  37  (^}  1  +  “7“  [“  37  (hlcos  v  > 


h  h  cos  v 

+  i—  (-7-2 - \ 

3z  '  s 


h  g  h  h  h  cos  X  3h 

“  7F  (h3cos  X  }  1  +  —  2 -  [77  “  37  (hlcos  v  >  1 

s 

(A-19) 


h^h2  g  h^^cos  X 

C7  =  ”7”  tcos  X  37  ( - e - 


2  2  2 

~  h0h0  h,h0h0sin  X  . 

)  _  L_  (_2_3)  1  +  2 _ 3_  (S_, 

'  3z  1  s  n  3  3  z  h„ 


hlh2h3  3  9h2 

+  - 2~  [cos  X  (iip008  X)  _  77“ ^ 


(A-20) 


3h, 


2  2 

hlh2h3Y  3  s  ,  h1h2h3COS  X  ,  -3i 

C8  3  3  z  h  “  2  3  z  (hiCOS|j)  9x  1 

s  3  s 


Vo  v2 
hlh3 
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,9  ,  9  ,  ,  hlh2  f  9  ,hlh3,  3 

[Jx  (h2C0S  X  }  -  97  (hicos  V  )  ]  +  -g—  [cos  v  57-  (-5—)  -  cos  A  ^ 


(■ 


h..h0cos  y  h,  h 


13 


■}  ]  +  "V  (h2cos  v  )  “  cos  y  g~ 
s 
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(A-21) 
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2 

h-^h^  g  h-jl^cos  X  g  h^hgCos  y  h-jl^cos  v  3h~ 

+  “s  [cos  y  97  <T - >  -  cos  X  37  ( - i - +  - 2 -  [9F 

s 

9  ^  ]_^2  9  3 

-  g*^  (h2COS  X )  ]  +  — 2“  [cos  X  ^  (h^cos  u  )  -  cos  y  gy  (h^cos  X  ) 


9  h. 


3  3  ^2  ^1^2  3 

-  cos  X  g^  (h2cos  v  )  +  cos  y  g^-1  +  — —  [57  U^cos  y  )  -  g^ 


h  h  h.h _cos  y  .  hnh0cos  X 

*  —  <57  ( — ; — >  - c°G  v  57  < — 5 — 1 1 


h  h  h  cos  X 
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( A-22 ) 


_  hlh2h3  f  .2,  9  ,19  /  s  .  .  9  ,19  ,s  u  A  .  .  -  3 

C10  2  th2h3sin  x  3 x  ^  s  9  x  h.  )]  h2h3Y  9y  {s  9x  (h..  hlh2S  9z 

S  XX 
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2 
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h?  9hn 
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(A-23) 
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,  _  hlh2h3  .  .  2,  3  .13  ,s  ..  ,  ,  3  .13  ,s  „ 

'll  — th2h3Sin  X  3x  {s  3y  (h  )}  +  hlh3T  3y  {s  3y  (hT)}  + 
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2.  y-Momentum  Equation 
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+  D10V  +  D11W  +  D12U] 


(A-26) 


where  (a,b,c),  S,  (a  ,$  ,y ) ,  and  V  V  are  given  by  (A-2),  (A-3),  (A-5)  and  (A— 6) 
respectively.  The  (i  =  1,  6)  and  the  (i  =  1,12)  coefficients  are  given 
by 
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,  .  h„h  h  y  hnh0h0cos  p  .  3h„ 
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D' - p - 57  (hj>  *  —  lcos  V  57  < - 5 - >  -  57  (— i-1 1 


4 


+  hlh2h3  ,9_ 


~r~  l37  in3 


(hocos  p  )  - 


3h  hh 

— i-i  +  —£—2 
3  z  s 


3  , h3hl . 


[cos  11  37  (~i“) 


h  h  cos  p  h  h  h  cos  p  3h 

57  - >1  *  ~Zr -  [5i  <V°S  »>  -  53T1 


(A- 36) 


h2h1SjLn  u  5h3  ^  h2h3hia  ,1  3 _  ,s_4  l_  9h2,  h2hly  8h3 

2  3y  2  s  3  z  lhJ  u2  9z  1  2  3x 

s  s  3  s 


h  h  eosv  h 


8 -  [S 

3  z 

h  „hji, cos  p 

3  2  1 

r— 

2 

3  z 

s 

3  h^cos 

X 

)  ] 

h  h  cos  p  h  9h 

(h  cos  p)  +  f-  - )  -  —  r-1] 

1  3  z  s  s  3  x 


hh  h  3h  h 

(hpcos  X )  -  —  (h,  cos  p  )  ]  -  — —  [—  -r-f-  -  —  —  (h„cos  X  ) 
3y  I  s  s  3  y  s  3  z  2 


h2hl 


,h2hlcos  u  ,  3  .Vl 

( - j - )  -  cos  »  n  ( - j 


3  z 

cos  X 


3  3  ^2^3  3 

-  cos  p  £-£•  (hjoos  v  )  ]  +  — - —  [cos  X 


3  x 


^2^1  9 

)  ]  +  — —  [cos  p  ^  (h^cos  X  ) 


cos  X  7^-  (h^cos  p ) ] 


(A-37) 


hlhlh  t  3  h„h 

D6  " - 3  5?  (iT>  +  — 

s  3 

,2 

21  r3  .  3  ..  .  , , 

+  — —  [g^-  (h-^cos  v)  -  —  (h^cos  X)J  + 


»  h0h  .  h0h„ 

r  o  ,2  J  ,  3  ,  2  J  ^xn 

[cos  V  (— — )  -  cos  \1  g-^r  ( — —  cos  X  )  ] 


h0h0cos  X  3  h0  - 
tL  j _  r _ z  ~ 


3  x  3  z  1 


(h, cos  p )  ] 


h2h  cos  X  h2  h  Sh  hh  hhhcoso  3h 

— 5 -  'i~  57  ,hicos  "  >  -  —  TT  -  57  > 1  * - I -  >5y 
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hh  h  h  h  cos  X  h 

(h,cos  A),  *4^  ,_i|_  ,Vos  a,  — )  -44 


3z  '"2l 


■^2^3 


8 


D7S  —  lcos"E  (—  008  11  >  -51  (— > 1  *  -J 


8  ,hih3 , ,  h2h3hlsln  v  )  s 

ax  h7 


hlh2h3 


ah. 


-l  3  w  J"L  o 

+  2 . .  [cos  p  JJ  (h^cos  p)  -  g-^2] 

s 

2  2 

h0hJi1a  ~  h-,h0h0cos  y  . 

r\  _  J-  o  .S  .  _L  <c  3  r  ®  /i_  \ 

d8 - 3  97  ( h7} - 5 -  [37  (h2C0S  v }  - 

s  Is 


9hl.  h2hl  r3 
3y  1  "  2  [ay 


a  hph„  „  h,h9  h0h..cos  v 

-  (h2cos  v  )  1  +  -g—  [cos  X  ^  (-^— )  -  cos  p  ^  ( - - - )  ] 


h2hi  ra  .  w“i 

+  — ~  I- —  (h„cos  X  )  -  cos  v  - — 

a  z 


a  x  s 
ahn 


2  lax  '  3 


hoh3  2  9h1  9hi  hph-4  9h,  hphJi,Y  , 

-  -V-  [hQ  sin  p  r-i  +  luce  r-*-  +  -M-  +  -M  1  | 


D, 


9  s2  3  K  ay  “2“  3z  ’  3x 


s3  ax 


l^h^cos  X  9h^  g 


hlh2 


[cos  V  r—  ( 


g  h^h^cos  p 


S  —  lal~-  a7  (h3cos  y  > 1  +  s  —  v  az  ■  s 


-  h„hpCos  v  h„hp  ,  ah, 

-  cos  “  a5  (^f - 11  +  -¥  'll  'V°6  *>  -  5^1 

h2h3  a  a 

+  — —  tcos  p  g^  (h-^cos  v)  -  cos  v  ~  (h-jcos  p  )  -  cos  p  |—  (hpcos 
s  ^ 

h_h„cos  p 


ah  h„h  h„h„cos  v 

+  cos  v  ^2.]  +  fl-  f  JL2 


s  ^8x  ^  s 


)  _cosx  - )] 


hlh2h2cos  y  a  a 

+  - 2 -  ^3y  ^h3cos  ** )  “  94  (h2cos  v)j 


D. 


10 


 hlh2h3 


-2  a  .13  .  s 


^  'Vl  sin‘,‘  5y  {t  iy  <4»  +  h3V  h  <7  57  <%» 

a  i  a  o  h0  h  h  cos  v 

+  h2h3Y  aT  {S  ay  (h^)}  1  +  i“  [a¥  {  I -  (a7  (h3cos  x  }  "  a7  {r 


(h^cos  p  )  ] 

(A- 38) 


(A-39) 

^cos  p ) 


(A-40) 


X  ) 


(A-41) 


•J^COS  V  )  )} 
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h^l 


cos  y 


3h"  9  hl  3h2  9 


*  37  -  <5?  (hloos  u  >  -  93r»  *  37  {—  (ai~  -  57  (h3oos  x  1 11  *  57 


2 

hph^cos  X  ^  r»  rv  h~  9 hp  ,* 

{ - - -  (g^  (hlCOS  v )  -  u  (h3cos  \ ) )}  +  97  {i2  (air  -  ay  (hicos  v  > » 


.  h,h0cos  u  » 

+  {-^Li —  -9 

3x  1  s 


3h, 


(ay  (h3cos  x)  "  a7 


=J  ] 


(A-42) 


Vi 

„  1  T3  ,L  ,  J  3  ,13  ,s  „  ,  u  ,  3  fn  ,s  „ 

D11  2  [h3hls±n  y  ay  (s  9z  (h  )}  h2hl  3z  {s  3z  (h  )} 

s  j  2 

1  1  8  r1  9  ,s  ,n  ,  h2  .a  h2hicos  v  8h3  a  a 

h2h3Y  37  {7  37  (h7)}  1  +  i“  l37  { - i -  <Tx  -  37  (hicos  y  } )}  +  37 


h^hn  cos  y  a  a  hn  .  3h„ 

s -  (3?  (hlcos  w  }  -  97  (h2cos  X  >  »  +  37  {~  (37  (h2COS  X  ]  ~  dT)} 


a  h/^h/^cos  X  ~  3ii«  r,  h 

+  37  {^~i -  (37  (hicos  y }  -  7T)}  +  37  {i“  (37  (h2cos  x  > 


h  h  cos  y  3h„ 

ly  (hicos  v  ) )}  +  fy  { - I - (ay1  -  h  (h2C0S  x  } )}  1 


(A-43) 


V\  Vn  V\ 

1  2  3.,,  .2  3  ,1  3  ,s  „  ,  3  ,13  ,s  „ 

Di2  =  — 2  (h3hisin  »  ay  {7  37  (h")}  +  h2hia  37  {7  37  (h7  } 

s  11 


3  ,13  ,  s 


+  h0h„y  £■—  {—  7—  (r—)}  ]  +  —  [3—  { 

2  3  3x  s  3x  n  s  3z 


h2  ,3  ,h2V°S  V  ,3 


3h. 


37  lh3oos  u  >  -  3i-)l 


g  h^h-^cos  y 


+  37  { 


3h, 


Igyi  -  |y  (h2cos  v  )  )}  +  |y  (|y  (h2cos  v  )  -  |y  (h3cos  y  )  )} 


.  h„h„cos  X  3h,  .  -  h„  .  3h, 

+  37{ - i -  (37“-  37  (h3cosu))}  +  37  {7"  (s7  (h2cosv)  -3 r)} 


h  h  cos  y 

+  37  {^7 -  (37  (h3cos  W  5  ‘  37  (h2cos  V  }  )}  ] 


(A-44) 


105 


3«  z-Momentum  Equation 


£  y_  W  3  W  13  uw  13vw  13  2  ~Z 

hx  3x  h2  3y  +  h"  77  +  h“  TT  +  hT  33T  +  h~  7i~  +  Tl™  +  T2(U  +  u 


(T  +  a)  uw  +  T  (V2*  v2)  +  TW2  +  (T  +  c)  w2  +  T  (UV  +  uv)  +  T  W 
^  4  4  5  6 


+  (T6  +  b)  vw  =  -  [sin2v  +  +  v  [V2W  +  En 


PS" 


h3  3  z  p  hx  3  x  '  “  h^  7y 1 


1  3  z 


3W  3W  9U  _  3U  ,  3U  3V  3V 

E2  3x  E3  3y  E4  3z  +  E5  3Y  +  E4av  +  E'7r7+ErtTr  + 


3  V 


J3  3y  ^4  3  3  *5  3x  e6  3y  +  E7  77  +  E8  77  +  E9  3?  +  Elo' 


>W 


+  E11U  +  E12V1 


(A-45) 


where  (a,b,c),  S,  (ct,g,y),  and  V2W  are  given  by  (A-2),  (A-3),  (A-5)  and  (A-6) 
respectively.  The  Tj  (1  =  1,  6)  and  the  E^  (i  =1,12)  coefficients  are  given 

by 


T  -  a  r  1  3  /■.  .  ,  1  3 

T1  ~  ^2  [hj^  37  (h2C0S  X  }  ~ 


h2hi 3y  (hlCOS  P 1  "  97  3y  + 


3h, 


_1 _ 3  ,, 

7^77  (h2C0S  v  }1 


(8  -  cos  u  sin2v  )  3hl  (sin2v  -  g 


hJi-jS 


3  z 


h 


3hlS< 


.  3h„ 
cos  p  )  _ 3 

!  3  x 


T  —  ®  [3  ,, 

2  77  ~2  [77  (h2cos  v }  - 


3h 


T3  = 


h2hlS 


sin  v  .3 


9y 

3h 


1,  sin2v_  f3hl  3 

2  Tz-  "  37  (h3cos  V  )  1 


hJ^S 


2  [37  (h3COS  X)  “  77 


h^S 


6 


3h 


2  3  ,,  .  . 

-  KT~  (xln  COS  V  )  ] 


rp  _  a  r3 

4  ■  7  .  c2  [7z  (h2cos  x }  - 

n^n^o 


3h, 


3y 


h^s2  3x  "9y  1 

3h, 


^ _ _  (•-  Z3  ^ _  /y.  \  i 

2  3  x  ~  77  (hlcos  V  )1 


h^h^S 


a  9h?  9hi  ,  2 

5  .  .  q2  3  x  "  C0S  V  3y  ] - , 

h^h^S  J  h  h  S 


3h 


12 

•  2 

sm  v  ,3 


3  1 


2  [cos  v  tt1  -  It  (h3cos  x)] 


+  .  ,  „2  ^3y  (h3coslJ  )  ~  3  z  (h2cos  v)) - p  [cos  v 


h3h2s 


h2hlS 


3  h  3  h 

— ^ _ L] 

3  x  3y 


(A— 46) 


(A- 47) 


(A-48) 


(A-49) 


(A- 50) 
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tz^lco 


cos  X 


3 


T6  = 


.  2  3h 
sm  v  f  3 


h3h2S 


2  3y 


3h2 

77 


3h 


2  3 


3h_  3  h 

a  r  .  3  _ £_■ 

o  lcos  x  9y  “  9z 


hlh2S 


[cos  X  g-J—  -  7jy  (h-^COS  u)J  + 


h^h-^S' 


2  lh  (hlcos  v  }  -  57  (h3C0S  X  ]  1 


(A- 51) 


2  2 

1  g  h3h!h2  2 

E1  =  "  iKJ  {3¥  (— ' S 


sin  v  )  +  yy  ( 


3  ,h3hlh26,  3  ."KV,,  *i3*1l*12s*n^v  3 

s — 1  *  w  (— s  »  *  —  s. 


~3  3  z 


h^gh^cos  y  9 


,gy  (h1cos  v)  - 


3h„  h..hJi,,cos  X 
_ 2,  12  3  .9 

3  x  J 


3h, 


2 -  [37  (h2cos  v  }  -  Tf 


h0h„  ,  h„h..cos  X  Lh  cos  y  h0h 

<C  f  a  ,  <£  1  x  O  ,  L 

*  8  ,C0S  v  37  1 - i - >  -  37  < - i 


-n  ,  s  .W03  “ , 

*  -j-  [COS  V  ^ - > 


.  h0hncos  X 

—  (-^= - )] 

3y  s  n 


(A- 52) 


2h^h„cos  X  .  h  h  h  cos  v  3h„ 

E2  =  ~2 -  lIy  (hlCOS  u  )  “  37  (h2cos  X  }  ]  +  - 2 -  (h2cos  X  }  "  37“] 

S  s 

hJi  g  3  h?h3  3 

+  -y2  [J-  (h^h-^cos  v)  -  cos  X  gy  (h^h-jcos  y  )]  -  — y  [h2  gy  (h1cos  y  ) 
s  s 

2  2 

g  h2hgCos  X  g  hJi-|h20  g  g 

+  h3  37  (hlcos  v  }  ]  +  -  2 -  37  (h2h3 C0S  X  }  +  3  37  (h7) 

s  s  ^ 


13  ,h3hlh26, 

s  3  z  s  ' 


(A- 53) 


2  2 

hArhuct 


3  3  z 

(h  ' 

s 

3 

vA 

h~h,cos  y 

ri_ 

2 

3  x 

s 

h_h^h_a 

13  ,31 

"  ) 

h0h. 


a  ^3^1  2  ^^3  g 

(h2cos  X  )  -  yy  (h1cos  y  )  ]  +  — y  [cos  y  yy  -  yy  (h2cos  X  )  ] 

s 


12  3 


3h, 


2 - [37  (hicos  y  ]  ~  tt1 


(A- 54) 
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2  2  2 

h„h,  h0sin  v  h„h„  .  h0h, cos  v  .  h„h, 

e4 - 3 - 87  V  ♦  —  loos  v  97  < — s - >  -  s7  1 1 

s  1 


hlh2h3  ,3 


3h2,  .  h3hl 


2 —  lyy  (h1cos  v)  -  yy— ]  +  — j [cos  v  ^7  ( 


3  ,hlh2, 


3y  v  s 

o 

9  ,hlh2eos\,  h3hlVos  v  ,9  „  ,  3hl, 

5?  <  5  >1  +  — 2 I37  (h2cos  v)  -  a^T1 


(A-55) 


hyijBln^  8hj  h^B  13  t  81,^  h^a  9hj 

5  -  -  s2  >*  ,2  '*H  S’  'h2  3x  -  s2  9y 


h  h  cos  X  h  h  h  cos  v  h  3h  h  h  h  cos  v 

-J-i -  [-2-  L-  fh  eo„  v\  4.  £_  I.J.  2 - )  _2 _ ii  123 

s  s  9  x  in2cos  v;  9z  1  s  1  s  9y  J - 2 

s 


3  3  hp  3  h..  ^  ^ 

[3l  (h3cos  M  5  "  37  (h2cos  v  >  ]  '  ~^T~  [7~  37"  ‘  7~  37  (h3C0S  u  }  ~  37 


h?h„cos  y  iuh?  3hp  huh 

(— — - )  ]  +  2  [cos  V  gTF-  ~  cos  v  a7  (h2cos  X  )  ]  +  y  1  [cos  y 


3_ 

3y 


,h3h2 


COS  V 

8  — )  -  COS  V  ( 


g  h^hpcos  y  h3h2 


)  ]  +  - —  [COS  V  gy  (h^cos  y  )  -  cos  y  ;y 


(h..  cos  v ) ] 


(A-56) 


h^hpa  g  h3h2  3  h3hi  9  ,h3h1cos  y 

E"  "  "5 - aT  +  “ —  tcos  X  yy  ( — — )  -  cos  v  —  (-* — 7 - 


"6  3  3x  vh_  '  s 

s  1 

h3hl  ,3 


3  x 


)  ] 


h^h-jCos  y  3^  g 


—  [gy  (h2cosX)  -  gy  (hyosy)]  +  - y -  lW  “  37  (h2C0S  u  >  1 


h  h  cosy  h  h  3h  h_h_  h,h0h.cos  v  3h_ 

*  -Jur-  h  )  -  t2  ^  -  h  (4^)1  -  -^4—  ■  1 


s  3y  3y  s 


3  z 


3  , ,  h3hl  ,h13  „  .  3  , 

-  37  (h3cos  y }  1  +  [F  37  (h3COS  U  ]  +  37  (" 


hJtycos  y  hn 


i — >  -  "r  h  (h2eos  v  i 1 


(A-57) 
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h0hn  ~  h0hncos  v 

•o  i  J-  r  d  ,  <d  1 

E7  =  T"  lcos  u  8?  < - 5 - 


)  - 


a_  ,h2hi% . 

9y  s  +  3 


h^-^sin  v 


—  (— ) 
9y 


hlh2h3  3  9hl 

+  - 2  tcos  v  37  (hlcos  V)  '  37" 

s 

h^lh26  3  ,s  ,  h1h2h3cos  v  (3 


e8  =  s3  ay  'h2 


(i~)  - 


5—  (h^cos  X  )  -  » 
9y  3  9  z 


(A-5 8) 


8h2,  h3h2  ,3 

(hn  cos  v  ) 


2  3  z 


^  h. 

-  9-33  (h3cos  X)]  +  -i— 


3  ,h2h3  3  h^008  x  h3h2 

[C0S  11  37  (_F"r  -  cos  v  37  (~^i - >  1  +  — 2" 


3 

[^7  ^hicos  V  )  “  cos  x 


(A-59) 


hJi  3h„ 

E„  =  -  —4^—  [Lsin  v  t —  +  hJ5 


2  ‘“l”"1  3z  3  3  x  h_  3y 


9h2  +  h3hla  9h2,  +  h3hlh2a  3  .  s  +  h3hlcos  y 


(— ) 
s3  9 y  ^ 


3  h 


2  3 


9x  3y 
.  .2 

h3hl  r3  , 

+  — —  [933  (h3cos  X  ) 
s 


h^tu  h.h  cos  v 

<L_  t-JL± _ 


2  1 

(hq  cos  v  )  ]  +  —  z  [cos  A  ( 

I  S  dX  S 


h  h  cos  X 
>  -  COS  U  37  (  s - 


)] 


9h2  hJ^  9  9 

9-3—]  +  - ■-—  [cos  v  9-33  (h2COS  X  )  -  cos  X  ^33  (l^cos  v  ) 


3h  h  h  h  h  cos  X 

■COSV7  (hlCOS  u  )  +  COS  X  gy-]  +  ~f~  [333  (-2 - - - 


h  h  cos  V  h  h  h  cos  v 

COS  ”  ty  (^4 - >  ]  +  "2 -  (l7  {hlCOS  V  > 


933  (h3cos  X  )  ] 


(A-60) 


_  hlh2Vl3  .  .23  ,13  , s  .  .  .  3  ,13  ,s  „  ,  ,  3  ,13 

E10  2  Ihlh2S1D  v  3  z  ^  s  9  z  (h„)}  hlh2S  3x  {s  3z  (h.)}  h3hl°  3y  {s  3z 

s  3  j> 


h  h  h  cos  X  h  h  cosy 

(h^)}]  s  [37{  s  (gy  (hlCos  y  )  -  933  (h2cos  X  )  )}  +  933  { - 3 -  (f^ 

2 

3h„  h?  3h„  h-h-cos  y 

<V°S  X  *  -  IF”  +  57  {i“  <3T  -  3?  <V°S  “  >  »  *  37  { - 5 -  <37  (h2oos  x  > 
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3y  '1 


3  hl  9N  3 

(h.,  cos  y  ) )}  +  —  {—  (5-^ - —  (h2cos  X  )  )}  +  |y  {  -  "s -  (jjy  (h1cos  y  ) 


3  rh2hlCOSX  ,3 


3y  s  '3y  3  z 


3h, 

3x 


-)}  1 


(A- 61) 


■n  _  12  3  r-i  ,  .2  3  r  1  3  f  s u  ,  ,  Q  3  r  1  3  .  s  u  ,  ,  3  r  1  3 

E11  2  Ihlh2sin  V  3z  {s  Dx  (L)}  h3h2B  3x  {s  3x  (h„)}  h3hl“  3y  {s  3x 

s  1  2 


h~  a  h~h?cos  X  3h,  » 

1  +  IT  <37  {  s -  <37-  -  Sx  (h2C0S  v  >  »  +  57  < 


3  ,hlh2COSV 


37  (h2cos  v 


3  .  .  „  3  :h2  ,3  ,  .  8hln  3  V1!00*3  11  3 

-  37  <h3cos  u))>  +  7T  {7~  (7T  (h3cos  p  ]  "  S7")}  +  aF  {  s - (37  (h2C0S  v  > 

2 

3h,  .  h,  .  .  „  h,h„cos  v  3h, 

"  9F)}  +  37  {T  (37  (h3COS  M  )  "  37  (h2cos  v  > )}  +  37  { - s -  (37* 


ay  (h^  cos  y  ) )}  ] 


( A- 62 ) 


hlh2h3  f,  ,  .23  fl  3  / s  vi  ,  ,  „  3  rl  3  ,s  „  .  ,  3 

E12  =  ~2  [hlh2Sin  v  37  {7  37  (h^)}  +  h3h2S  37  {7  37  (hj)}  +  h3hla  37 

,1  3  ,s_u  ,  h3  r3  fh3h2COS  X  3  9h2n  3  h!h2C0S  v  3h2 

s  3y  (h2)}  ]  s  3x  ^  s  (3y  (hlcos  v  }  ~  3x  )}  +  3x  {  s  (3z 

2 

3  a  h?  3  a  3  h^h-,  cos  y  3  h? 

-  ay  (h3cos  X  ) )}  +  ay  {—  (37  (h3cos  X)  -  37  (hlCos  V  )  )}  +  ay  {-J— i -  (37- 

2 

r.  .  h1  h  h.cos  v 

-  ay  ( Vos  V  ) )}  +  ^y  {_  (_  (h3cos  X  )  -  37^  +37  { - I -  (37  (hlcos  v  } 


gy  (h3COS  X  ))}  ] 


(A-63) 
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C.  Turbulence  Equations. 


8  rVt 


U  3k  V  3k  W  3k  _  1  .3  ,  t  3k  .  „  3k  A  _  3ku  A  0  .  „ 
^  57  hj  3J  +  h~  37  -  ¥  [77  {7^  {A  37  +  H  37  +  G  37)}  +  37  *777 


/«  9k  +  R  3k  .  -  3ku  3  /t  lr  3k  .  _  3k  .  _  3kn  .  *  /A  ... 

(H  37  +  B  37  +  F  37)}  +  77  {7T  (G  37  +  F  37  +  c  “)}  1  +  ^  (A"64) 


3  z 


U_  §£  *  I_  i£  A  W_  3£  _  1  r3_  /t  3e  Je  Jen  3  rVt  3e 
h^  3 x  h2  3y  h^  3z  s  ^3x  ^£s  ^A  3x  +  H  3y  +  G  3z^  +  3y  B  3x 


+  B  fy  +  F  fz)}  +!7{7T(Gl7+Fl7+cl7»)  +Gei7G-Ge2F 

e 


(A-65) 

where  s  and  the  coefficients  (A,B,C,F,G,H)  are  given  by  (IV-13)  and  (IV-14) 
respectively.  The  turbulence  generation  term  is  defined  by 


G 


*  v2(4*  4 


)  +  4(e 


2 

12 


)] 


( A-66 ) 


where  e 


ij 


is  the  rate -of -strain  tensor 


eij  - 1 17i  *  ’xT) 


(A-67) 


T 

In  (A-67)  VV  is  the  deformation-rate  tensor  e^j  and  V V  its  transpose,  i.e., 
T 

VV  =  e^.  The  components  of  are  defined  by: 


'll 


=  e 


11 


2  (h1AA11+  h1HA21+  h1GA^1) 


12 


s 

1  , 

2  e12+  e21 


'21 


1  ^1  ^2 

2  [—  (AA12+  HA22+  GA32)  +  (HAn  +  BA^  +  FA31)J 

s  s  J 


(A-68) 


(A-69) 
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(A- 70) 


el3  2  (el3+  e3l}  e3l 

1  B1  ^3 

=  2  [—  (AA13  +  HA23  +  GA33)  +  -f  (GAn+  FA, +  CA  )  ] 
S  s  ^ 

^2 

e22  "  ~2  (HA12+  BA22+  FA32} 


£23  2  (e23+  e32} 

1  ,h2 


'32 


h. 


2  I—  (HA13+  BA23+  FA33)  +  —  (GA12+  FA22+  CA^) 


P  =  P 

33  33 

h„ 


— 2"  (GA^3+  FA23+  CA33) 


where 


(A- 71) 

(A- 72) 

(A-73) 


An  =  V  Uan*  Va2i*  Wa3i 

Ai2*  Dbn*  V  ">31 

V  Ucn+  V=21*  *2  W=31 

A21=  V  Ds12+  Va22+  Wa32 
a22=  ub12'>  V2+  vb22+  *b32 

V  H°12*  VC22*  V  "C32 

A31=  V  Ual3+  Va23+  "a33 

V  “13*  V  Vb23+  *33 
A33=  Dc13*  VV  V  *C33 


and  the  notation 


,,  ...  ,  ,3V  3 V  3W  , 

(Bi>  Vi>  Wi>  -  <3ir.  it:- 

ill 


(A- 74) 


has  been  used.  The  coefficients  (a^-,  b-y,  cij- )  in  (A-74)  are  the  components  of 
the  vectors  representing  the  derivatives  of  the  unit  vectors  e.,  in  the  (x,y,z) 
coordinate  directions,  i.e. 
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(A-75) 


ei,j"  8Xj  (ei) 


aijei+  bije2+  °ije3 


where 


hi  ■  (air  hi(1,/hi 

bn=  bii'hi 

cll= 

al2  '  {&12~  hl,2)/hl 

b12"  h^^l 

C12  = 

a13  =  (al3"  hl,3)/hl 

b13=  b13/hl 

C13= 

a21  =  *21^2 

b21=  (b21-  h2, 

l)/h2 

C2l" 

a22  =  a22^2 

b22=  (b22‘  h2. 

2)/h2 

C22  = 

a23  =  a23/h2 

b23=  (b23”  h2, 

3  )/h2 

°23= 

a3l  =  a31/h3 

b31=  b31//h3 

C31= 

a32  =  a32/h3 

b32=  h^2/h3 

c32= 

a33  =  a33/h3 

b33=  b^/h3 

c33= 

'11 

,  i 

'12 


/h 


C13/hl 

C21/h2 

C22/h2 

C23/h2 
(c3l“  * 

(c'  -  i 

^32 

(c'  -  V 

1  33 


3,1 


(A- 76) 
)/h0 


3.2) /h3 

3.3) /h3 


and  the  notation 


has  been  used  with 


8x,  (V  ~  hi,j 


(a 


11'  U1P  ull' 


=  D 


-1 


hlhl,l 

v  *  i  “  ^1^1,2  ^ 

*',1  “  hlhl,3 


(A- 77) 


(a 


22'  22'  °22 


)  = 


_1  h2h2'2 
D  ]  X  ^2  ”  h2h2,3 
V ’2  "  h2h2,l 


(A-78) 


fo  1  h  1  c  •  ) 

a33’  33'  33 


,-l 


h3h3,3 


"  D  U  u',3  h3h3,l 
x  ^3  -  h3h3,2 


(A-79) 
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and 


(a12'  h12>  c12 )  D 


-1 


(a 


13>  13’  13 


Cio)  -  D 


-1 


(a23'  h23’  c23)  "  D 


-1 


hlhl,2 

¥2,1 

\  ^',2~  v',3+  X)l)  J 

1  hlhl,3 

^  (X  1 1—  ]i  1 2 +  v  1 3) 
h3h3,l 


T  (v'o-  X  '  +  y  '  ) 


2  "-,3 
^2h2,3 
h3h3,2 


,2' 


-1 


D 


hlh2h3S 


v  1  =  h^t^cos  v 


(A-80) 


(A-81) 


(A -82) 


hph^sin  A  h^h^y 

hlh2S 

h2h3y 

2 

h^h^sin  u 

hxh2a 

( A—  83) 

h2h36 

h-jh^ot 

h^h2sin^v 

A  ' 

=  h^h^cos  A 

U  ' 

=  h^hyjos  m 

(A-84) 

Note  that  the  coefficients  (a^,  b^,  c^)  are  symmetric,  i.e.,  (a^J,  b^,  c|j) 
(aji>  bji'  Gji}  for  1  *  J* 
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APPENDIX  II:  Partially-Parabolic  Equations 


In  Section  IV-A  the  procedures  for  obtaining  the  partially -parabolic  form  of 
the  governing  equations  have  been  discussed.  The  resulting  equations  are 
provided  in  this  Appendix.  The  definitions  of  the  various  coefficients  appearing 
in  the  equations  are  given  in  Appendix  I. 


x-Momentum  Equation 


U_  9U_  V_  9U  W_  9U  2  2  1_  9uv  1_  9uw 

hl  9X  +  h2  9y  +  h3  9z  +  *  +  h2  9y  +  h3  92 


[  s  in^A 


2_2l>  2^  in 

\  Y  h2  3y 


+  8 


1_ 

h. 


9  z 


s 


[— 

9y 


(a 

a22  9y 


9 _ 

9y 


ft  ) 


3_ 

az 


(a 


ajj 

32  9  y 


9 

3  z 


,  9JJ.  . 
a33  9  z  1 

( B— 1 ) 


y -Momentum  Equation 


U  3V  .  V  3V  .  W  9V  .  1  9v  .  1  9vw  ^  3  U  ..  ^  3U  A  ,e  ^  — 

D?  sy  -  ”  %  Jz  *  <E1+C)  w 


+  8.2'  +  Sj  ( U2  +  u2  I  +  ( £_;  *  b)  v2  +  S5  (UW  +  uw)  +  SJJV  +  (S^+a)  uv 


1  r.2  19jD  1  3i  A  1  9_p,  v  ,3 

-  - -  — ^  +  a  — ^  +  y  —  — ^]  +  —  [- — 


[sin  u 


9  V>  3  .  9  V. 

a0o  — )  +  r~:  (a00  r— )  + 


pS* 


h0  9y  h„  3 z  1  h,  3xJ  s  l3y  VQ22  9y'  3y  va23  9z 


9 _  .  9V  9_  9V 

9 z  a32  9y  +  3 z  la33  9z' 


(B-2) 


z -Momentum  Equation 


U  9W  V  3W  W  3W  1  9w  1  9vw  „  9  U  „  3  U  m  Tnir 
hx  3x  +  h2  3y  +  h^  9z  +  h^  3z  +  h2  3y  “  V  E6  9y  vE4  8z  Tl™ 
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+  T0  (U2  +  u2)  +  (T  +  a)  uw  +  T„v2  +  (T .  +  c)  w2  +  TC(UV  +  uv) 
2  1  3  4  5 


+  (T^+  b)  vw 


1  r  *  2  1  IE.  a  q  1  3P  A  1  ap,  ^  v  r9  r 

“2  Isin  v  hT  at  +  B  TT  +  «  h7  a£]  +  s  [3t  <a22  3y} 

po  5  1 


+  L_  (a  il_ 
3y  la23  9  z 


3  .  3  W.  3  ,  3W,  , 
3 z  (a32  3y  3 z  (a33  3z}  1 


(B-3) 


k-Equation 


vB  vF  \>F  v  C  ' 

U_3k  +  V_3k  +  W_3k=I  ,3_  ,_b_  3k  _t_  3k.  3_  (_L_  3UL  +  _L_  3k,  i 

h  3x  h_  3y  h  3z  s  3y  a  s  3y  a  s  3z  3z  a  s  3y  a  s  3z 

J.  X  K.  K.  K. 


+  G  —  e 


(B-4) 


e  -  Equation 


U_  3e_  V_  3e_  £_  3e_  =  3.  f3_  V  tB  3e_  ^  3e_.  3_  VtF  3e_  Vl  3e_. 

h  3x  +  h  3y  +  h  3z  s  3y  a  s  3y  +  a  s  3z  +  3z  a  s  3y  +  a  s  3z 
123  e  e  e  e 


(B-5) 


The  Reynolds  stress  terms  required  in  equations  ( B— 1 )  -  (B-3)  are  related 
to  k-e  through  the  isotropic  eddy  viscosity  concept: 


v.v. 
i  J 


-  2v ,  e  .  . 
t  ij 


*  3k  (hih,1Rij) 


( B— 6) 


where  v 

(A-73), 

tensor 


2 

.  =  C  k  /e  is  the  eddy  viscosity,  e  is  the  rate-of -strain  tensor  (A-68)- 
t  u  ij 

h^  are  the  metric  coefficients  (IV-10),  and  g.jj  is  the  inverse  metric 
(IV-14.1 ) .  Note  that  equations  (B-l)-(B-3)  were  derived  under  the  assump- 
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tion  that  v.v.  ~  0(e).  Thus,  the  order-of -magnitude  of  v  e..  ~  0(e);  however, 

1  J  X  lj 

the  order-of -magnitude  of  v  is  difficult  to  assign  and  should  be  0(e  )  <  v 

t  —  t 

2  -1 

0  (e  ).  For  this  reason  terms  of  0(e  )  and  0(1)  have  been  retained  in  the 

partially -parabolic  form  of  e...  Note  that  e..  is  also  required  for  the 

t  J  i  j 

evaluation  of  the  turblence  generation  term  G  (A-66) .  The  terms  (A-74)  simplify 
to 


11 


^2  = 

^11 

A13  = 

Ucll 

A21  = 

U2  + 

A22- 

^12 

A23  = 

UC12 

A31  = 

U3  + 

A32 

Ub13 

A33  = 

OA 

i — 1 

O 

t=> 

°i +  Uan 


12 


(B-7) 


and  the  components  of  e . .  (A-68)  -  (A-7 8)  become 

ij 


^1 

E11  =  2  ^A(ui+  Uaii>  +  h(U2+  Uai2)  +  G(U3+  Uai3)} 

s  ^  A 


(B-8) 


1  hl 

e12  =  t  [~t  (AUbu+  H((Jb12+  V2)  +  G(Ub13+  V3)}  + 
s  ^  * 


h2 

+  -|{H(U1+  Uan)  +  B(U2+  Ua12)  +  F(U3+  Ua13)} 
s  ^  ^ 


(B— 9) 


1  hl 

e13  =  2  1  2  {AUcn+  H(Uc12+  W,)  +  G(Uc13+  W  ) 

S 


h3 

+  (GtU^  Uan)  +  F(U2+  Ua12)  +  C(U3+  Ua13)} 
s 


(B-10) 
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(B-ll) 


e22 

E23 

e33 


2 

=  —  {HUbn+  B(Ub12+  V2)  +  F(Ub1^+  V^)} 

1  rh2 

-  2  f—  {HUcn+  B(Uc12+  W2)  +  F(Uc13+  W^)  +  {GUbn+  F(Ub12 

s 

+  C(Ub13+  V  )}  ] 

h3 

-  —  {GUcu+  F(Uc12+  W2)  +  C(Uc13+  W3)} 


(B-12) 

(B-13) 
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APPENDIX  III:  Equations  in  Discretized  Form 


This  appendix  provides  the  details  of  the  finite-difference  procedures 
used  to  put  the  partially -parabolic  Reynolds  and  Turbulence  equations,  as  well 
as  the  pressure-correction  and  pressure  equations,  into  discretized  form.  A 
written  description  has  already  been  provided  in  Section  IV.  B  and  IV. C.  and 
will  not  be  repeated  here.  As  shown  below,  various  terms  in  the  equations 
have  been  grouped  for  convenience.  All  differences  are  labeled  as  to  forward 
(FD),  backward  (BD),  or  central  (CD).  Lastly,  the  terms  are  collected  and  the 
coefficients  in  equations  (IV-18)  -  (IV-22)  and  (IV-29)  -  (IV- 30)  are 
defined.  The  following  definitions  are  used  for  the  distances  between  nodes 
(see  figure  3) : 


(Ax  )4  = 

+  m,n 

[(X4+1 

-  x4  )2  + 

(/+i 

i 

ru 

i 

m,n 

m,n 

m,n 

m,n 

(Ax  )4  = 

[  (X4  - 

X4"1)2  + 

(Y4  - 

Y4"1)2 

-  m,n 

m,n 

m,n 

m,n 

m,n 

m,n  m,n 


m,n  m,n 


(Ay  )4  =  [(Y4  -  Y4  )2  +  (Z4  .  -  Z4  )2]1/2 

v  ^  +  m,n  m+l,n  m,n  m+l,n  m,n 

(Ay  )4  =  [(Y4  -  Y4  )2  +  (Z4  -  Z4  )2]1/2 

*  -  m,n  m,n  m-l,n  m,n  m-l,n 

(Az  )4  =  t  (Y4  ,-  Y4  )2  +  (Z4  -  Z4  )2]1/2 

+  m,n  m,n+l  m,n  m,n+l  m,n 

/*  _  r +  ftp £  _  £  *2,1/2 

-  m,n  m,n  jn,n-l  m,n  m,n-l 

(Ay  )„  (Ay  ) 


--  m,n  -  m,n 


/A  (A  \4-l 

Az  =  (Az  )  ^ 

—  m,n  -  m,n 
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where  (£,m,n)  are  the  indices  in  the  (x,y,z)  coordinate  directions  and  (X,Y,Z) 
are  the  Cartesian  coordinates.  Note  that  in  the  present  staggered  grid 
arrangement  the  streamwise  velocity  component  U  is  located  at  the  grid  node 
(Z ,m,n)  whereas  the  transverse  velocity  component  V  is  located  at 

(Z -1/2,  m-1/2,  n),  the  girthwise  velocity  component  W  is  located  at 

A 

(£-1/2,  m,  n-1/2 ) ,  and  the  remaining  variables  (p,p,k,e)  are  locate  at 
(£-1/2,  m,  n).  All  the  required  geometrical  quantities  are  first  evaluated  at 
the  grid  nodes  (£,m,n)  using  a  backward  difference  for  x-derivatives  and 
central  differences  for  y-  and  z-derivatives.  Subsequently,  when  the  discre¬ 
tized  equations  are  formed,  all  geometric  quantities  are  evaluated  at  the 
location  of  the  variable  under  consideration  by  taking  the  appropriate  average 
of  the  neighboring  values  which  is  the  reason  that  the  notation  £  ±  1/2,  m  ± 
1/2,  n  ±  1/2  has  been  introduced  below. 

A.  Reynolds  Equations 

1.  i-Momentum  Equation 


a) 


if- 


!_  3U  _  jrun _  (u£  _  u*-l)  =  xa  (U*  _  u*-l)  (BD) 

h  9x  ,  ,£  v  m,n  m,n'  1  m,n  m,n 

1 


(C-l) 


b) 

V9 


L  -L-{(V1)  0»m+  -  l*y -1)  Cl,n>  -  Ibltl,n-  *>2^ 

V 

+  xb^l/"  _ 
i  m-l,n 

v*  -  I  'C.l  *  Cl.n  +  Cl.n1 

Hec=  Re  V*(4y++  sy.)‘>n/2 


b23y  Ayv 


where 


(C-2) 
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Re  -  l/(vt+  1/RJ, 

=  -  1  and  Ayv  =  2(AyJ*  n  if  ReQ  >  2  (BD) 

<j>y  =  0  and  Ayv  =  (Ay+  +  AyJ*  ^  if  | Re^ |  <  2  (CD) 
<f>y  =  1  and  Ayv  =  2(Ay+)^n  if  Rec  <  -  2  (FD) 


c) 

wau _ 

h„3z  Az 

3  w 


W 


((V  11  <„*1-  *  'V1'  “!n,n-l>  -  ”1 


(C-3) 


where 

* 

W 

Re 


(/  +  /  +  l/+1  +  /+1  )4 

m,n  m,n+l  m,n  m,n+l 

W*  Re  (Az  +  Az  )£  /2 
+  -  m,n 


=  -  1  and  Azw  =  2(AzJm  n  if  Rec  >  2  (BD) 

4>z  *  0  and  Az^  =  (Az++  Az  )*  n  if  |Rec|  <  2  (CD) 
<f>  =  1  and  Az  =  2(Az,)^ 


w 


+  m,n 


if  Rec  <  -2  (FD) 


d) 


R4u  +  2  h23y  (vte12"  3  k12)  +  2  h^z  {vte13”  3  k13)  "  R4r2+  2{h^3y  (vte12' 
ik12}  +  h^a"z  (vte13"  I  kl35  *  h^Ty  (vt  h^y)} 


R 


u2  +  -  8 


2(i 


o  o )  ^ 


(V^E,,-  4  ^o) 


4  Vy  *  12  3  12  h3az  *  13  3  13 

-  2  ls4(<V^7E12,i*l,n 


(vts12"  3  k12*m-l,n  1  *  4z  {(vtE13"  3  k13lm,n-H"  <'ltE13“  3  k13lm,n-lt  1 


xd  (CD) 


( C— 4) 
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2 


Vy 


(v 


au 


t  h28y  (Ay  +  Ay  )~  (Ay  T 

+  -  m,n  +  m,n 


(  r 

t  —kV.1/?!  n  (Dl  _  rf  J 

£  , .  ,£  'm+l,n  m,n 


(v.) 


-  m~1/2’n  (ui  -  ui  1  )}  =  -  xdp  U*  +  xd~  0*  -  xd .  U*  , 

j£  m,n  m~l,n  2  m+l,n  3  ra,n  4  m-l,n 


(CD) 

(C-5) 


where 


-  3  U 

e12  e12  ”  2hp3y 

V 

4ys  ■  (4yt  *  4y-C'n/2 
4ze  =  (4V  4z-tn/2 


e) 


2 

1_  .sin  X  3p  y  3p  g  3p 
P  1  g2  1^3  x  +  g2  hp3y  +  g2  lu3z' 


2  p£+1-  p* 

1_  ,  ,sin_X_.£  pm,n~  pm,n 

p  1  c2  V.n  Ax 
S  ’  x 


(1_)£ 


1  .  £+1  £  £+1  £  S  £  1 

Ay  (pm+l,n  +  Pm+l,n“  Pm-l,n  “  Pm-l,n^  +  '  2'm,n  Az 
X  o 


.  £+1  £ 

pm,n+l+  Pm,n+1 


£+1 

p 


1 


£ 


)} 


xe1  (FD),  (CD),  (CD) 


( C— 6) 


where 


4*x  -  ' 

4yx  =  (4y+*  4y->l,n 

Az  =  (Az  +  Az  )^ 
x  +  -  m,n 


f) 
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v  3  3  U.  3  ,  3  U  3  .  3U  3_  3U„ 

?  {3y  (a22  dj]  +  3y  (a23  3^}  3¥  (a32  3y)  3z  (a33  3z)} 


(I)  +  (II)  +  (HI)  +  (IV) 


v  3  ,  3U. 

(I)"  o 


s  ay  '“22  ay  (h  h  s)*  (4y  ,* 

1  3  m,n  +  m,n 


(h2a22)m+l/2,n  (u*  _  /  }  _  (h2a22)m-l/2,n 


m+l,n  m,n 


(Ay  ) 


(if  -  (f  i  J>  - 1 

m,n  m-l,n  ,  + 


{A*-  +  Ay+}m,n 


_  =  xfn  if  -  xf„  if  +  xf.if  ..  (CD) 

£  1  m+l,n  2  m,n  3  m-l,n 


(C-7) 


/  T  T  \  _  v_  3 _  /  3  U'  _  r 

( II ) -  g  (ar:  a  rr '  I 


3y  la23  3z'  "  lh1h3S(iy_-  %n 


m+l,n 


_  (  h3  ajl_)£  (if  -  if  .)}  =  xf,  (CD),  (CD) 

^Az  +  Az  ,  v  m-l,n+l  m-l,n-l  4 

+  -  m-l,n 


(C-8) 


3  (J  ^  ^  d  a  o 

(III)=;slz  (a32  3y}  =  {h1h„s(Az  +Az  f  {  (Ay_+Ay+)m,n+l  (  m+l,n+l"  m-l,n+l) 

X  *  —  in,n  “ 


(a2  32  f  (u* 

'Ay  +  Ay^  ,  m+l,n-l 

J-  t,+ra,n-l 


-  if  t  n  -,)}  =  xf.  (CD),  (CD) 


(C-9) 


,T.n_  V  3  ,  3U»  _  r  2v 

(  IV)  -  s^rr(a^'3^r7>)  t 


az  a33  az1  ‘h^S  (4z_*  4^)'%n  (40+)«^ 


t  (h3  a33)m,nH/2  J.  , 

J  I  £  '  m,n+l  m,n 


(h~a 


l  33  m,n-1^2  (ufc  _  u*  )}  =  xf  (CD),  (CD) 
£  '  m,n  m,n-l  6 


(C-10) 


(Az  ) 

—  Ill,  IX 

Collecting  all  terms  (C-l)  -  (C-10): 


«1<„  -  +  “itl,.-  IVi,n*  *  XC1  *  xdl 

"  xd2  4+1  ,n  +  xd3Ura,n"  xd4Um-l,n  *  xel  "  xf  l4i+l,n  +  xf2Um,n  xf3Um-l,n 


-  xf,  -  xf  -  xf  =  0 
4  5  o 
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or 


(xb3  -  xd4  -  xf3)  U*_M  +  (xax  -  xb2  +  xd3  +  xf2)  U^n  +  (Xbx  -  xd2  -  xf^ 
U4  =  xa^U4'1 


-  . n  =  xan u  -  xc,  -  xd..  +  xf,  +  xfc  +  xf,.  -  xe1 
m+l,n  1  m,n  1  1  4  5  6  1 


Finally,  the  coefficients  in  equation  (IV-18)  are: 


a,  = 

xb  .  -  xd , 

-  xf 

1 

3  4 

= 

xa,  -  xb„ 

+  xd 

2 

1  2 

= 

xb.,  -  xd„ 

-  xf 

3 

1  2 

Su  = 

xa1U4-^'  - 

xcn 

1  m,n 

1 

Pu  = 

-  xej 

(C 


2.  y-Momentum  Equation 


a) 


U9V 

h^3x 


U 


(Ax  )  ,  n 


“T75 —  (V4m  -  V4"1)  =  ya^V4  -  V4-*)  (BD) 
-!/<;  m,n  m,n  1  m,n  m,n 


(C 


where 


U*  =  (U4  +  u£'1+  U4  n  +  ir4-* 

m,n  m,n  m-l,n  m-l,n 


)/4. 


b) 


-11) 


-12) 
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V* 

K^y  ■  iff  f'V11  il.nVv  -  y»l  Cl.n-^.n*  ^-l.n 

(C— 13) 


where 


=  -  1  and  Ay_  =  2(Ay)*“^  if  Re  >2  (BD) 

V  m-l/2,n  c 

=  0  and  Ayv  =  (Ay++  Ay_)^'^  n  if  |R®c |  1  2  (CD) 


$  = 


.4-1/2 


=  1  and  Ayy  =  2(Ay+)m_1/2  R  if  Rec  <  -  2  (FD) 


Re  =  Ee  /  (Ay  +  Ay  Y'1^  /2. 

c  m,n  -  +  m-l/2,n 


c) 


wav  =  w 

h09z  Az 
3  w 


{  (<b  +  1)/  .-2 |  /  +  (<j>  -1)  /  ,}  =  yc, 

L'vz  m,n+l  yz  m,n  vvz  '  m,n-l  J  1 


(G— 14) 


where 

$ 


W 


4  -1/2 

-  1  and  Az  =  2 (Az  )  ,,,  „  ,  /0  if  Re„  >  2  (BD) 

w  -  m-l/2,n-l/2  c 


+  z  ■  0  “d  4zw  -  (4a+*  4z-Cli/2,„-l/2  lf  l»»el  ‘  2  (CD) 


z 

* 


1  and  Az  =  2 (Az  )^ „  if  Be  <  -  2  (FD) 
w  +  m-1/2 ,n-l/2  c 

(/+/.+/,  +  W2,  -  _  )/4. 

m,n  ra,n+l  m-l,n  m-l,n-l 


Re 


*  2 
ReW  (A  z++  Az  _)  *  /2. 


d) 


2{ 


(v^e.J  + 


(v.  e„0-  4  kn  )} 


h  9y  '  t  22'  h  3z  '  t  23  3  23 


2  9k22  _  0!  a 

3  h2ay  " 


(Vj 


h^y  t  2h2ay 


av  J  +  8  (vte22) 


(V  .  8/ 


1  -  u  2  9k22 
k00)}  -  -o 


h  3z  '  t  23  3  23"  3  h23y 
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28  (vte23~  3  V  2  8  ^22  28  (vte22} 

h^z  ~  3  h29y  +  hgSy 


Az  {  {vte23  7  k23)m,n+l+  (vte23~ 


(k„^  -  5* 


1  ir  — 1  _  1  —  ,  lr  >1-1  i  2  22m, n  22m-l,n 

7  k23  m,n+l"  (  t  23"  3  k23)m,n-l~  (vte23"  3  k23)m,n-l}  "  3  Ay  + 


,£-l 


Ay  ^Vte22^m,n+  ^Vte22^m,n  ^vte22^m-l,n  '  vt°22  'm-l,n 


/  - /-I  , 

-  (v  e„„ )  J 


=yd1(CD),  (BD),  (BD) 
(C-15) 


3  ,  3  V  , 

(V+  ^-V")  =  ~ 


^-i/2 
^  o  m,n 


h  3y  t  h  3y'  .  ,£-1/2  1  £-1/2 

2  2  (Ay  +  Ay  )  (Ay  )  ' 

+  m-l/2,n  +  m-l/2,n 


(V*  -  /  ) 

m+l,n  m,rr 


(v 


,£-1/2 


(477^f-  Vtl,n»  -  -  ^3<n-  ^/m-l,n<0DI 

1  y-Vl/2,n 


(C-16) 


where 

Az 


,  .  ,£-1/2 

=  (Az  +  Az  ) 


Ayc  = 


'22 


/ .  ,£-1/2 
l4y->m,n 


1  3  v 


'22 


2h2  8y 


e) 


_  1  (_I _ -3-P  ,  Sin  y  3p  a _  3^  1  ,^£-1/2  1 

n'9  Jv  O  o-,r  O  ‘  «  11  iL  I/O  ^  lu 


_ _ _  _  _  /  £  +1 .  £  +1  £  -1 

p  's^  8x  "  S2h  8y  ^  S2h  8z/  ~  p  1  lg2^m-l/2,n  Ax^  pm,n  pm-l,n_  pm,n 

12  3 

2 

_  A-1  %  +  fsln  u  ,£-1/2  1 _  ,  £  A  *  .  ,ot_,£-l/2  1  .  I  £ 

pm-l,n  2  ;m-l/2,n  Ay  lpm,n  pro-l,rr  (  2}m~l/2,n  Az  (pm,n+l  pm-l,n+ 


"  Pm,n-l“  Pm-l,n-l 


)} 


ye1  (CD),  (BD),  (CD) 


(C-17) 
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where 


AXy  =  {Ax++  (Ax_+  Ax++)/2} 

a  /a  ,4-1/2 
Ay„r  =  (Ay  )  _ 


m,n 


.  ,A  .  ,4-1/2 

Az  =  (Az  +  Az  ) 

y  +  -  m-l/2,n 


f ) 


v  ,3  ,  9V,  3  .  3  V,  3  .  3  V,  3  ,  3V„ 

s  {3y  a22  3y  +  3y  (a23  3z}  +  3z  (a32  3y5  +  3z  (a33  3z)}  =  (I) 


+  (II)  +  (III)  +  (IV) 


/  T  v  _  v  9 {  3  V\ 

(I)  s  3y  a22  3y 


2v 


{ 


,h  o  -1/2 

<tZ  nun 


(h„ano) 


4-1/2 


rv  x.  o/»  .  ,4-1/2  L  ,A  ,4-1/2 

O^hyStAjr.*  4y+  m-l/2,n  ,ly+,m-l/2,n 


(V*  ,  -  /  ) 

m+l,n  m,n 


:  %Vl'D  ^,n-  VLl,„»  *  yfl  ^l,n-  yf2  yf3  /m-l,n(CD) 

(C— IS) 


,4y-)m-l/2,n 


,TT,  _  v  3  ,  3 V, 

(II)  s  3y  (a23  3  z 

,,  ,4-1/2 

h0a00  ' 

3  23  m-l,n 


fVl  Q  ,4-1/2 

-i  (Ji0a0  0 ) 

_ v _ 1_  .  3  23  m.n _ 

.  o  ,4-1/2  Ay  1  /A  A  ,4-1/2 
(hlh3S)n.-l/2,n  a  (4z.*Az-  Vl/2,n 

-  (V*  ,  V* 


(V*  _  v^  ) 

vm+l,n+l  m+l,n-l 


,  r/o  .  i  —  r  .  .)}  =  yf .  (CD),  (CD) 

(Az  +Az  )A  m  1,n  1  m  ^,n_^  ^ 

+  -  m-l/2,n 


(C-19) 


where 


4ya  =  {4y-  *  I  ,4y*+4y-,}m‘n/2 


,TTT,  _  v  3  .  3 V, 

(III)  s  3 z  (a32  3y 


,h  n  ,4-1/2 

_ v _  ,  1  2  32  V-l/2.n+l  fv4 

,  ,,4-1/2  1  Ay,  lvm+l,n+l 

{ h, h0S(Az  +  Az  )}  _  h  ’ 


12 


J-  m-l/2,n 
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,  .£- 1/2 

/  ,  _  Zi  a32)m-l/2,n-l  /  _  / 

vm-l,n+l;  Ay  m+l,n+l  vm-l,n-la 

G 


=  yf5(CD),  (CD) 

( C— 20) 


where 


4y„  -  {4y_*  |  <4y*+  4y-»m"n»l 


4yc  '  {4y-  *  7  l4y.*  4y-»l'n-l 


(XV)  =  —  —  (a  — ) 
1  ;  s  3 z  la33  3 z 


2v 


£-1/2 

(h,h?S)  '  (Az  +A  z  )  ,  „ 


(h  a 

_  r  3  33Vl/2.n>l/2 

£-1/2  1  £-1/2 

4 z+  m-l/2,n 


(/  /  ) 

m,n+l  m,n 


,,  .£-1/2 

n3&33  m-l/2,n-l/2  /  /  n 

.£-1/2  Vm,n”  Vn-l" 

{Az-: )m-l/2,n 


=  yf6  (CD),  (CD) 


(C-21) 


g) 


{ (S-^+c)  ( — 2v^e 23^  ^  ^23^  +  ^2 ^  3  k33'  ^ S^(U  — 2v ^x  +^^X1^ 


+  (S^+b)  (^  k22~  2vte22^  +  ^5^  ”  ^vte13+  3  k13^  +  ^UV  “  ^^S6+a^vte12  ” 


1  t:  u£-l/2 
3  k12)}m-l/2,n  7gl 


( C— 22 ) 


h) 


vft7  H  *  uD8  H  -  vD7Z/2>n  tU^.n-  Cl,n  *  C  '  CJ.n’/2 


+  vD8£_1/2  ((/  _1  +  U +  lt“3 


m-l/2,n  m,n+1  m»n+1  ®~l,n+l  m-1, 


n+1 
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(C-23) 


-  U*  /-1  ,  -  /  ,  ,  -  U*  }  n)/8  =  yhn  (BD),  (CD) 

m,n-l  m,n-l  m-l,n-l  m-l,n-l  '  J  1  ’ 


Collecting  all  terms  (C-12)  -  (C-23): 


ya,  (/  -  b)  +  yb,  ,  -  yb„V^  +  yb0V^  n  +  yc,  +  yd, 

1  m,n  m,n  J  1  m+l,n  J  2  m,n  0  3  m-l,n  J  1  0  1 


-  yd2  Cl,n  +  yd3^,n  -  ydX-l,n  *  yel  '  ^iCl.n  *  **£,n 


-  yf3  ’i-l.n  -  yf4  -  yf5  -  yf6  +  ygI  -  yhl  =  0 


or 


(yb3  -  yd4  -  yf3)  /_1>n  *■  (yaj  -  yb2  *  yd.,  *  yf2)  />n 

*  lyV  yd2‘  yfl>  Cl,n"  yaX’J-  yol-  ydl*  yV  yV  yf6~  yRr  yV  yhl 


Finally,  the  coefficients  in  equation  ( IV— 19 )  are: 


-  yd^  -  yf 


b2  =  ya-j^  -  yb2  +  yd 3  -  yf2 


(( 


b3  =  ybl  *  yd2  -  yfl 


/-1 


Sv  =  yal  V  yCl"  ydl+  yf4+  yf5+  yf6"  ygl+  yhl 


Pv  =  -  ye 


-24) 
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3.  z-Momentum  Equation 


a) 


U9W 

h  &x  ' 

1  (Ax  ) 


\-i/r - 

-'m,n-l/2 


za,  (wj,  -  wj,  i)  (BD) 
1  m,n  m,n 


where 


U*  =  ((/  +  l/  1+  l/  1+  (/■  ,  )/4. 

m,n  m,n  m,n  m,n-l 


b) 


(C-25) 


h29y  Ay^  "ry 


m+l,n 


2d»  w  +  (+  -l)  vr  ,  } 

vy  m,n  VHy  '  'm-l,nJ 


=  zbn 


(C-26) 


where 


£  -1/2 

*y  •  -  1  and  4y„=  2<4y .)m>n.1/2  If  Rec  >  2  (BD) 

*y  =  0  and  ayv  =  (4y+*  if  |ReJ  <  2  (CD) 

*y  =  1  and  4yy  =  2(Ayt)‘^2/2  if  Reo  <  -  2  (FD) 

-  *>  ^n-1/2/2- 

V>  "  <n  +  Cl,n  *  <,n+l+  Cl,n*l>"- 


c ) 


«i,n*l  -  2*Z<,n  *  (♦,-!) 

X  TUT  7  7  7 


W 


=  zc-.v/'  _  -  zc-l/'  +  zcJ//'  i 

1  m,n+l  2  m,n  3  m,n-l 


(C-27) 


where 


£  -1/2 

h  ■  -  1  “d  4zw  =  2(l2-)n,,n-l/2  lf  Reo  >  2  (BD> 


*z  =  0  and  4zw  =  <4V  ^z-^m~n-l/2  lf  lHeJ  <  2  <CD) 
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£-1/2 


0,z  =  !  and  Azw  =  2(AZ  +  )ffl>n_1/2  if  Rec  <  -2 

Re  =  Re  /  (Az  +  Aza)£_1/?/0 
C  m,n  -  +  m,n-l/2 


(FD) 


d) 


^  _  2  8 

(V  +  eOO“  ^OO  )  +  v:  (^  +  £Qo) 


2  9  *33 


h_9y  '  t  23  3  23  h.  3z  'vt°33'  3  h.9z 


g  ^  _ 

—  2  U  ^  „  (^4^0*)“  knn  )  ■*" 


h23y  t  23"  3  23 


h^3z  (vte33)  “  h^3z  vvt  h^3z;  “  3  h^z 


_iw_}  _  2  d_hi 


p  (—5 (v  e  A  >  )  +  — (v  E  n  _ 2A  _  _A f  (v  e  —  k 

^^y  'vte23'  3  23'  h^z  1  t  33 ^  3h^9  z  Ayg  Uvt  23  3  23V+l,n 

+  (v  e  _  A  _(v  e  _  i  V  -(ve  -  —  k  }  -  — 

1  t  23  3  23;rn+l,n  lvte23  3  *23'm-l,n  lvt  23  3  *23  Vl,tf  3 


,r-  £  r-  £ 

(k„„ _ -  k 


n£  -1 


%£  -1 


^>n  ,3>,nrl  +  {  —  ,t  (v^00)l  _  (v+e00)r" 

£-1  Az  t  33  m,n  —  — 

(Az  ). 


t  33  m,n  '  t  33  m,n-l  '  t  33  m,n-r 


zd1  (CD),  (BD),  (BD) 


(C-28) 


where 


Ay  =  (Ay_+  Ay+) 


£-1/2 

m,n 


Az.  =  (Az  +  Az^)^  /2 


33 


=  e 


+  m,n 
1  9W 


33  2h3  9z 


h^9z 


(v 


9W 


t  h^9z 


{ 


(vY^2 

u  in  j  n 


r>  ,  ,£-1/2  x  IA  ,£-1/2 

{ Az  +  Az  }  (A z  )  ,  n 

+  -  m  „  -1/0  +  m,n-l/2 
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,  >£-1/2 

— —  m;n~l _  lyf'  _  yfi  n 

(Az  )l~1/21/0  m’n  ^ 

-  m,n-l/2 


-  zd2<,,n*l*  *d3<i,iT  zd4<,n-l  (CD)  (C- 291 


e) 


.  2 


1,3  9p  a3p  sin  v  3p 
"  p  l02  9 x  +  n2^  +  02. 

1 


j  =  _  1  ,(B_j4-l/2  1 {  4+1  +  4+1  _  4  -1 

S“ hn  S^y  Sfch39z  p  S2  m'n~1/2  Axz  m,n  m,n-l  pm,n 


£-1  ,  a_  £-1/2  1 _  (  £  £  £  £ 

^m,n-l  ^2  m,n-l/2  Ay  ^m+l,n  -^m+l,n-l  ^m-l,n 


.2  £-1/2 

*  h-  (rm,n-  *  zel  <CD)-  <CD»-  (BD) 


( C— 30) 


where 

Ax 


\  {Ax  +  j  (Ax  +  Ax  )}* 

^  ^  +  m,n 


Ay„  = 


- .  ^  .  >4-1/2 

(Ay++  Ay-)m,n-l/2 


Az  =  (Az) 


4-1/2 


f ) 


v  /L_  fo  Mi  tL  9W>  A  9  ,  9 W  9  ,  9W„ 
s  {9y  a22  9y  3y  (a23  9z*  3z  (a32  9y}  +  9z  (a33  Fz)} 


=  (I)  +  (II)  +  (III)  +  (IV) 


,-n  _  v  9  ,  9W, 

(I)  s  9y  a22  9y 


(h~a„„ ) 


4-1/2 


_ 2 _  r  2  22  m+1/2  ,n-l/2  .<Z  JL  . 

lhh3Si^-K'Xl/2  W^Cn-?/2  "tl,B  m’n 


.  .4-1/2 

m2  a22  m-1/2  ,n-l/2  ...JL  ..JL 

-(~>l/2  m,n~  Vl,n»  =  zfl(CD) 

,y-  m,n-l/2 


(C-31) 


/tt\  _  v  9  /  9W N 

(II)  -  a23  Tz  ' 


,,  .4-1/2 

_ I _  f  3a23  m+l,n-l/2  * 

<V>3S  (4y_+  4Z»  m*1'n*1 
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(hoaoo) 


£  -1/2 


_  y/'  )  _  _ 3  33  m-l,n-l/2  ^JL  n  _  _f  (CD)  (CD) 

m+l,n-l  A2l  m-l,n+l  Vl,n-l,j  zf2tCD)'  (CD) 


(C— 32) 


where 


i2a  =  !Az_*j  (Az/- 


4%  {Az_+  ^  tAz++  Az — ^jn-*l^n 


/  T  T  T  ^  ^  ^  ^  W  % 

(III)  s  az  (a32  ay 


(h  a  )*  -1/2 
^n2a22  m.n+1/2 


(W^  - 

'  "m4-1  1^4.1 


,,  v»  ca  »A-l/2  .  x£-l/2  v  m+l,n+l~  m-l,n+l' 

.  ,*-1/2 

-  (<n,n-r  ■  zf3  (CD)-  (0DI  (C-33) 

(4y-*  ^'.n-l'n-l/a 


where 


Az  =  {Az  +  \  (Az  +  Az 

c  -  ^  +  —  m,n 


(IV)  =  1 —  /a  iLK) 
UVJ  a z  la33  a z 


(h_a00) 


£-1/2 


,£-1/2 

_ 2 _ (ria33.V,n  (w£  Jl  , 

/,  v,  o  a  ,£-1/2  1  .  ,£-1/2  "m,n+l  wm,n 

(h  h_S  Az  '  /0  Az  )  '  ’  ’ 

12  a  m,n-l/2  +  m,n-l/2 


1*1  n-  ^  n  i  ))  =  zf,/  n+1-  zf,/  +  zf,/  ,  (CD) 
*£-1/2  m,n  m,n-l  4  ni,n+l  5  ro,n  6  m,n-l 


{Az-)m,n-l/2 


(C-34) 


where 


4  1  ,4  4  ,£-1/2 
Aza  =  ^  (Az/  4z-)ja>n— 1/2 


g) 


'T1W  *  VM  hi-  2Vll>  *  (V«H-2vtA13,§i13)  ♦  T3(f  k22-  2v 22 ) 

+  ( T^+  c)  (j  k-o-  2vte;33)  +  T5(UV  -  2vte12  +  |  k12 )  +  (T6+  b) 


(- 


2vte23+ 


—  V 

3  23  ;m,n-l/2 


zg-. 


( C— 35) 
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h) 


au  „  9U  „  £-1/2  1  .1-1  1 

vE6  3y  +  vE4  3z  vE6m,n-l/2  ra+l,n+  m+l,n  +  m+l,n-l 

+  U2-1  -  U2  -  U2-1  -  U2  .  -  tf'}  )/8 

m+l,n-l  m-l,n  m-l,n  m-l,n-l  m-l,n-l 

+  ve/-1/?/0  (U4  +  -  U4  „  U t~l  i  )/2  =  zh,  (CD),  (BD) 

4m,n-l/2  m,n  m,n  m,n-l  m,n-l  1  ’ 


(C-36) 


Collecting  all  the  terms  (C-25)  -  (C-36) : 


za,  (W2  -  W2  "*")  +  zb,  +  zc^W2  ,-  zc-W2  +  zcJN2  n 
1  m,n  m,n  1  1  m,n+l  2  m,n  3  m,n-l 

+  zd,  -  zdJIlf2  .,+  zdJK2  -  zd.W2  ,  +  ze,-  zf,  -  zf,,-  zf., 

1  2  m,n+l  3  m,n  4  m,n-l  11  23 

-  zf.W2  ,+  zf.W2  -  zf,/  .+  zg,  -  zh,  =  0 

4  m,n+l  5  m,n  6  m,n-l  tol  1 


or 


( zc _—  zd,  -  zf,)/  ,  +  (za,  -  zc„  +  zd„  +  zf-)  / 

3  4  6  m,n-l  1  2  3  5  ra,n 

+  (zc,  -  zd.  -  zf.)  W2  =  za.  W2-'*'  -  zb.  -  zd.  +  zf,  +  zf,,  +  zf, 
12  4  m,n+l  1  m,n  1112  3 


-  zg1  -  ze][  -  zhx 


Finally,  the  coefficients  in  equation  (IV-20)  are: 


c„  =  za,  -  zc„  +  zd„  +  zf 
2  12  3  5 

c^  =  zc1  -  zd2  -  zf  (C-37) 

Pw  =  -  ze^ 

Sw  =  za.  W2-1  -  zb.  -  zd,  +  zf,  +  zf-  +  zf„  -  zg,  +  zh, 

1  m,n  1112  3  11 
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B.  Turbulence-Model  Equations 


F  stands  for  k  or  s 


a) 


=  S!L  (F2,  -  F*_1)  =  Fa,  (F*  -  F*"1)  (BD) 

h^3x  Ax  ra,n  m,n  1  m,n  m,n 


(C— 38) 


where 


U*  =  +  U*_1)/Ax 

m,n  m,n 


1  i 
Ax  -  ^  {Ax  +  Ax  }  m  ^ 

2  -  —  m,n 


b) 


V  3F  =  Jh. 


h23y  Ay 


-  (/  .  -  F4  ,  )  =  Fb,  (F*-,  -  F*  ,  )  (CD) 

m+l,n  m-l,n  1  m+l,n  m-l,n 


(C-39) 


where 


,  s2-l/2 

Ay  =  (Ay++  Ay_)mjn 


c ) 


w  —  = 

9z 


W 


Az 


w 


{($  +1) 
z 


F*  ,  -  2<f> 
m,n+l  z 


F*  + 
m,n 


<V1} 


/  ,}  = 

m,n-l 


Fc, 


where 


4> 

4> 


z 


z 


2. -1/2 

=  -  1  and  Azm  =  2 (Az  )*  r  if  Re„  >  2  (BD) 
w  -  m,n  c 

=  0  and  Az  =  (Az  +  Az  )*_y2  if  |Rel  <  2  (CD) 
w  +  -  m,n  1  c 1 


(C-40) 
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£  -1/2 

*z  -  1  and  Azw  =  2(Az+)m>n  if  Rec  < 


-  2  (FD) 


d) 


Re  =  ReAzW 
c 


W*  =  (/  +  /  -|  )/2 

m,n  m,n+l 

.  ,4-1/2 

Az  =  (Az  +  Az  ) 

+  -  m,n 


1  ,3 


9F 


9F 


9  F, 


9  F 


7  %  tCv  a22  &  +  5y  <Cv  a23  +  ^Cva32  ^  +  97  <Cva33  37 


)} 


=  (I)  +  (II)  +  (III)  +  (IV) 


C  =  vi/a,  for  k 
v  t  k 

C  =  v  /a  for  e 
v  t  £ 


, _ ,  1  9  .  _  9F < 

(I)  =  s  9y  (Cv  a22  9y 


1  2  r  2  v  22  m+l/2,n 

r  4  -1/2  ,  ,4-1/2  1  ,  ,4-1/2 

{hlh3S}m,n  (A^+Ay-}m,n  (A^)m,n 


(h  C  a 

(F^  -  F^  )  -  2  v  22  m-l/2,n  ,^4  _  p4  )}  =  Fd  F^  -  Fd  F^ 

m+l,n  m,n  (Ay  ,4-1/2  1  m,n  Vl,n"  ralVl,n  ra2  m,n 


+  Fd-F8,  .  (CD) 
3  m-l,n 


(C-41) 


/  tt  \  1  3 t  n  3  F  x 

(II)  “  s  9y  (Cva23  3 z 


- 1 - ( (h3a?A)t~1/^  pi  , 

{h^Uy.^X^  4“-*’*  m+l,n  m*1'n*1 


h0a_0C  4-1/2 


-  ( 


■p—-r~V-)  .  (F^  .  .-  F4  ,  -.  )}  =  Fd .  (CD),  (CD) 

Az  +Az+  m-l,n  m-l,n+l  m-l,n-l  4 


(C-42) 


.Ttt\  _13  , n  3F , 

(III)  "  7  9z  (Cva32  9y) 


h  a  C  4  -1/2 

I  (  £  ^  v ) 

'  A  IAtt  fy 


{ h, h,S(Az  ♦  Azt)>*-1/2  4y--4y*  m>n+1 

12  -  +  m,n 
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( C— 43) 


( TT^ 

v  m+l,n+l"  m-l,n+l ' 
=  Fd5  (CD),  (CD) 


h  a  C  1-1/2 

-  (2  32  V )  (re 


-  / 


Ay  +Ay+  m,n-l  m+l,n-l  m-l,n-l 


)} 


(IV)  =  —  - _  (C  a  — )  = 

uv'  s  dz  1  va33  9z 


{(- 


h  a,  C  £-1/2 

(Fi.vrfi.) 


-) 


r  i«  v  o  /  *  a  \  ^  1  /  2  A  z  ..  /  ^  m .  ri  m  >  n 

{h1h2S(Az_+Az+)}m  n'  +  m,n+l/2  ,  ' 


_  (h3^33Cv)<l-1/2 


»,  '  (F*  -  F*  .  )}  =  Fd,  (CD) 

Az  ,  ,,  m,n  m,n-l  o 

m,n-l/2  *  ’ 


( C— 44) 


G  -  e  =  Sk 


Cls  l«  -  C2e  t  ■  & 


u  ,2 


,  C  k  £ 
£  ,  v 


i)  For  k,  G-e  =  G-  -iL  k  =  (G)*  -•  (-=— )  k*  »  Fe,-  Fe-k* 

7  v,  m,n  v,  m.n  1  2  m,n 

t  9  t  m,n  } 


(C-45) 


ii)  For  e,  C,  f:  G  -  C~  -  -r—  =  (C,  f  G)£  -  (C,  fr/  /  =  Fe,-  Fe0  / 

*  lc  k  2e  k  le  k  m,n  2e  k  m,n  m,n  1  2  m,n 


Collecting  all  terms  (C-38)  -  (C-46): 


Fa,  (/  -  F£_1)  +  Fb,  (F£  ,  -  F*  ,  )  +  Fc,  -  Fd,  /  , 

1  m,n  m,n  1  m+l,n  m-l,n  1  1  m+l,n 

+  Fd,/  -  Fd.,/  n  -  Fd,-  Fd,-  Fd,-  Fe,  +  Fe,/  =  0 
2  m,n  3  m-l,n  4  5  6  1  2  m,n 


or 


■<Fbl*  Fd3>  'i-l.n*  (fV  Fd2+  Fe2»  (Fbr  Fdl>  4*l,n 


(C-46) 


-  Fc  +  Fd,  +  Fdc  +  Fd,  +  Fe, 
1  4  5  6  1 
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Finally,  the  coefficients  in  equations  IV-21  and  (IV-22)  are: 
(d,e)1  =  -  Fb1  -  Fd^ 


(d,e)2  =  Fa^  Fd2  +  Fe2 

(C-47) 

(d,e)  =  Fbx-  Fd1 


SF  =  -  Fc.  +  Fd,  +  Fd_  +  Fd,  +  Fe. 

1  4  5  6  1 

C.  Pressure -Correct ion  and  Pressure  Equations  The  discretized 

the  momentum  equations  ( IV— 18) — ( IV-20)  can  be  put  in  the  form 

form  of 

/  =  0*  +  —  Pu£ 

m,n  m,n  m,n 

(C-48) 

/  =  V*1  +  “  Pv£ 

m,n  m,n  m,n 

( C— 49 ) 

/  =  w*1  +  —  Pw£ 

m,n  m,n  m,n 

(C-50) 

where  the  pseudovelocities  are  defined  as 


t,n  *  -  ^  (alCl,n*  a3Cl,n  '  <„» 

(C-51) 

t  =  -  TT  1  +  b-V*.,  -  Sv£  ) 

m,n  1  m-l,n  3  m+l,n  m,n 

(C-52) 

W*  =  -  —  (c,/  ,+  e/  Sw£  ) 

m,n  c2  1  m,n-l  3  m,n+l  m,n 

(C— 53) 

and  the  pressure-gradient  term  can  be  expanded  to  yield 


i;  <n  *  a7  Pm,n»  +  a8  <C,n*  4*1, n‘  V^rT 

/  4+l  A  4+1  4  , 

*  a9  Pra,n+1+  PB,ntl-  P»,n-1  '  Pm.n-11  (FD)-  ,CD)>  (CD) 

(C-54) 

1^4  ,,4+14+1  4-1  4-1  .  ,  .  4  4  . 

b2  ^nqn  “  b7^pm,n  pm-l,n~  pm,n~  pm-l,n^  +  b8^pm,n“  pm-l,n^ 

+  b9  {pm,n+l+  Pm-l,n+l"  Pm,n-l“  Pm-l,n-l)  (CD)’  (BD)>  (CD) 

(C-55) 
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(C—  56) 


1  rj*  _  ,  £+1  £+1  £-1  1-1  ,  £  £ 

c  ^m,!!  c7  %,n  pm,n-l"  pm,n  “  pm,n-l^  +  c8  *pm+l,n+  pm+l,n-l 


-4-l,n-pLl(n-l>  +  09'Pm,n-4,n-l>  (0D)'  (CD)-  ™> 

in  which  the  (a^,  b-^,  c^,  i  =  7,8,9)  coefficients  are  defined  by 


a7  a8  a9 


b?  bg  b9 


c7  c8  c9 


•  21 
sin  X 

Y 

8 

2  > 

S  4xxa2 

2  * 

S  Ayxa2 

S  Azxa2 

1 

sin  y 

a 

sHV 

8 

s24yyb2’ 

ot 

S24,yb2 

sin  v 

2  » 
S  Ax  c0 
z  2 

2  » 
S  Ayzc2 

2 

S  Az  c~ 
z  2 

(C-57) 


and  evaluated  at  the  grid  node  (£,m,n)  . 

Substituting  (C-48MC-50)  and  (IV-27)  into  (IV-26)  results  in  the 
desired  equation  for  pressure: 

1  ,  „  ,  l  ,~J l  ,  J l  ,  £+1  £  .  ,  %£  ,  £+1  .  £ 


(Ax  )2  I  (^2^3® )ra,n{ **m,n~  ^a7^m,n  ^pm,n_  pm,n^  "  ^a8^m,n^pm+l,n+  pm+l,n 
-  m,n 

-  Pm-ln-  4-1, n>  -  <Vm,n  ‘4Jn+l*  4,n*r  C-!'  4,n-l» 

-  [h  h  S)1'1  {^-  (a7)^  (p2  -  p^)  -  (aft)^  (p2  ,  +  v~]  . 

2  3  ro,n  m,n  7  m,n  m,n  m,n  8  m,n  m+l,n  m+l,n 

£  £-1  ,  ,  .£-1  ,  £  ^  £-1  £  £-1  , 

Pm-l,n-  Pm-l,n  a9  m,n  Pm,n+1  Pm,n+1-  Pm,n-1~  Pm,n-1 

+  ~  ,£-1/2  1  (hlh3S)m+l/2,n{^m+l,n~  (b8)m+l/2,n(pm+l,n_pl,n)  “  (b7)m+l/2,n 

1  y-  m+l/2,n 

-  £+1  £+1  £-1  £-1  ,  ,.  ,£  ,  £  £  £  £  ,, 

Pm+l,n  Pm,n-  Pm,n~  Pm+l,n  ~  °9  m+l/2,n  Pm,n+1  Pm+l,n+l-  pm,n-l”  pm+l,n-l" 

-  (h-h-S)2  w_  {V2  -  (bQ)2  . /0  (p2  -  p2  .  )  -  (bp.)2  1/0 

1  3  m-l/2,n  m,n  8  m-l/2,n  ^m,n  ±m-l,n  7  m-l/2,n 

.  £+1  £+1  £-1  £-1  ,  ..  ,£  ,  £  £  £  £  , 
Pm,n+  Pm-l,n~  pm,n-  pm-l,n  -  b9  m-l/2,n  pm,n+l+  pm-l,n+l-  pm,n-l~  pm-l,n-l  ^ 


(Az-)m,n+l/2 


-1/2  [  (blh2S)m,n+l{4,n+l"  lt!9)m,n+l/2(piii,n+r  ^n1  "  ^^n+l^ 
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^m+l,n+l  ^m+l,n  pm-l,n+l  pm-l,n^  -  ^c7^m,n+l/2^pm,n+  pm,n+l  pm,n  pm,n+l^ 

(hlh2S)m,n-l(*m,n'  <09lm,n-l/2lpm,n~  pm,n-l 1  "  (o8 )m,n-l/2 (Pm<-l,n  +  Pm*l,n-1 

'  pLl,n-  '  <c7>L,n-l/2<<n*  C-l'  C"  C-l»  '  =  0  «>-*> 

Through  the  use  of  the  following  definitions 

CX1  =  (^2^3^)m,n^(Ax-  ^n^n’  CX2  =  (h2h3S)m,n/(Ax-)m,n 

Cyl  =  (hlh3S)m+l/2,n/  {^ -] m+1/2 ,n>  cy2  =  (hlh3S)l-l/2,n/(Ay-)l+l/2,n 


.£-1/2 


CZ1  =  (^^S)m,n+l/2/(Aa- *111,11+1/2  *  Cz2  "  (hlh2S)m,n-l/2/(Az-)ni,n+l/2 

m  =  ex,!/  -  cx^l/  b+  cy,V^  n  -  cy0V^  +  cz,W^  . -  czJl/’ 

1  m,n  2  m,n  J1  m+l,n  J2  m,n  1  m,n+l  2  m,n 


alx  (a7  ^nqn'  aly  ^S^n’  alz  "  (a9]m,n 


-L 

7  m,n>  a2y 

=  (a8 

\*-4.  /  »X/  "*“_L 

m,n’  a2z  a9  m,n 

)l 

7  m+l^^n’ 

a3y  = 

(b8)m+l/2,n> 

ao  = 

3z 

(h  )£ 

v  9  Vl/2,n 

)*• 

7  m-1/2,11’ 

a4y  " 

^8^-172,  n’ 

a4z  = 

(b9}m-l/2,n 

)l 

7  m,n+l/2' 

a5y 

,  .£ 

°8  m,n+l/2’ 

a5z 

.  .£ 

°9  m,n+l/2 

)A 

7  m,n-l/2> 

a6y  = 

a6z  = 

f  / 

°9  m,n-l/2 

the  pressure  equation  can  be  put  in  the  form 

flCn  *  P2pro,n+  f3pm‘n  +  Vn*l,n*  f5pl-l,n*  f6pl,n+l*  f7pi,n-l  =  SP  ‘°-59) 
where 
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f,  =  -  a.  cx  -  a  cy  +  a.  cy  +  a,  cz  -  ac  cz. 

1  lx  1  3x1  4x  2  6x  2  5x1 

f .  =  a,  cx  +  a^  cx  +  a.  cy  +  a,  cy  +  ar  cz  +  a,  cz, 

2  lx  1  2x  2  3y  1  4y  2  5z  1  6z  2 

f  =  -  a  cx  +  a  cy  -  a  cy  +  a  cz  -  a  cz 

3  2x  2  3x1  4x  2  5x1  6x  2 

f  .  =  -  an  cx  +  a„  cx,—  a,  cy,-  ac  cz  +  a,  cz_ 

4  ly  1  2y  2  3y  1  5y  1  6y  2 

f5  ■  alycxr  a2ycx2-  a4ycy2+  a5y°V  a6yoz2 
f6  '  -alzoxl+  a2zcx2”  a3zcyl+  a4zoy2-  a5zozl 
f7  "  alzcxr  a2zcx2+  a3zoyr  a4zcy2~  a6zoz2 


and 


£+1 

(alzcxl+ 

£+1 

3xcyl’ 

■^m+l,n+ 

a5xral> 

^m,n+l 

£+1 

,aiycxl+ 

a4xcy2’ 

£+1 

6xcz2> 

Pjn,n-1- 

•^m-l,n 

(C-60) 


*  *a3z°yl+  a5yozl’  pm+l,n-H  <a3zoyl*  a6ycz2*  pm*l,n-l 

-  (a4zoy2*  a5yczl>  pm-l,n+l+  (a4zcy2*  a6yoz2>  pLl,n-l 

I  \  A  “I  /  i  A— 1 

-  (a„  cx„+  a0  cy, )  p  ..  -  (a,  cx0+  a,  cz, )  p  ,, 

2y  2  3x  ^1  *m+l,n  2z  2  5x  1  ^m,n+l 

,  ^  >  £-1  ,  ,  £ -1 
+  (a,  cx„+  a.  cz0  p  .+  a,  cx,+  a,  cy,)  p  , 

2z  2  bx  2  ^m,n-l  2y  2  4x  J2  ^m-l,n 


(C-61) 
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